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Abstract 


A  new  method  for  orbit  prediction,  which  is  as  accurate  as  numerical  methods  and  as  fast 
as  analytical  methods,  in  terms  of  computational  time,  is  desirable.  This  paper  presents 
Kolmogorov  Arnol’d  Moser  (KAM)  torus  orbit  prediction  using  Simplified  General  Per¬ 
turbations  4  (SGP4)  and  Two-Line  Element  (TLE)  data.  First,  a  periodic  orbit  and  its 
Floquet  solution  are  calculated.  After  that,  perturbations,  which  are  on  the  order  of  10  5 
and  smaller,  are  added  to  the  periodic  orbit  plus  Floquet  solution.  Then,  the  low  eccentric¬ 
ity  KAM  torus  is  least  squares  fitted  to  the  SGP4  and  TLE  data.  The  performance  of  the 
theory  is  presented  in  various  ways.  The  new  method  is  approximately  five  times  more  ac¬ 
curate  for  the  best  fits  and  three  times  more  accurate  for  mean  fits  comparing  to  SGP4  and 
TLE.  History  of  TLEs  and  KAM  torus  theory  can  be  used  to  make  accurate  orbit  predic¬ 
tions,  which  is  conceptually  similar  to  extrapolation.  In  addition,  the  new  method  may  rival 
numerical  methods  and  it  can  be  used  for  collision  avoidance  calculations,  and  formation 
flight  applications.  However,  high  eccentricity,  polar  and  critical  inclination,  air  drag,  and 
resonance  problems  should  be  addressed. 
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KAM  TORUS  ORBIT  PREDICTION  FROM  TWO  LINE  ELEMENT  SETS 


I.  Introduction 

Artificial  satellites  are  indispensable  parts  of  the  daily  lives  of  people.  Both  civilian  and 
military  world  depend  on  satellites  orbiting  Earth  at  very  high  speeds.  There  are  many  ap¬ 
plications  of  satellites  such  as  communication,  navigation.  Earth  observation,  and  weather. 
The  position  of  a  satellite  should  be  known  within  a  certain  limit  to  establish  communi¬ 
cation  between  ground  stations  and  the  satellite.  The  accurate  positions  of  satellites  and 
debris  are  important  for  collision  avoidance  purposes.  Orbit  determination  methods  yield 
the  position  of  a  satellite  at  a  given  time  in  the  future. 

1.1  Motivation  and  Background 

Bright  moving  objects  in  the  sky  have  intrigued  mankind  since  ancient  history  began. 
Ancient  people  worshiped  heavenly  bodies  as  Gods.  They  developed  timekeeping  systems 
to  facilitate  their  daily  lives.  They  used  stars  to  navigate  in  deserts  and  seas.  They  made 
observations  to  figure  out  the  motion  of  celestial  bodies.  It  was  Johannes  Kepler  who 
first  developed  the  three  laws  of  planetary  motion  in  1619.  Kepler’s  empirical  findings, 
based  on  Tycho  Brahe’s  planetary  observations,  were  explained  by  Isaac  Newton  in  his 
book  “Principia”  in  1687.  Isaac  Newton  used  his  three  laws  of  motion  to  derive  Kepler’s 
laws  of  planetary  motion.  Kepler’s  equations  serve  as  a  complete  solution  to  the  two- 
body  problem  (2BP),  which  determines  the  mutual  gravitational  interaction  of  two  bodies 
under  the  assumption  that  there  is  no  third  body  affecting  them.  The  solution  to  the  2BP 
assumes  that  there  is  no  perturbing  force  affecting  the  two  body  motion.  The  formulation 
of  the  three-body  problem  (3BP)  was  provided  by  Joseph  Louis  Lagrange  in  1772.  The 
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3BP  is  a  special  case  of  the  n-body  problem  (NBP).  It  is  the  problem  of  calculating  mutual 
gravitational  interaction  of  three  bodies. 

At  the  end  of  the  18th  century,  remarkable  advances  occured  in  perturbation  theory  due 
to  the  logarithms  invented  by  John  Napier  in  1614  and  the  efforts  in  modeling  the  Earth’s 
gravitational  field.  Perturbations  can  be  defined  as  the  divergence  from  the  normal  motion. 
The  motion  of  a  body  under  the  influence  of  small  accelerations  diverges  from  the  two 
body  motion.  There  are  three  different  methods  to  find  perturbations  to  the  two  body  mo¬ 
tion:  special  perturbation  techniques,  general  perturbation  techniques,  and  semianalytical 
techniques.  Special  perturbation  techniques  rely  on  numerical  integration  techiques  and  re¬ 
quire  high  computational  power.  Numerical  integration  methods  provide  orbit  predictions 
with  a  level  of  accuracy  in  meters.  General  perturbation  techniques  make  use  of  analyt¬ 
ical  solutions  to  the  perturbation  problem  and  do  not  require  much  computational  power. 
However,  predictions  made  by  general  perturbation  techniques  are  not  as  accurate  as  that  of 
numerical  methods.  Semianalytical  techniques  combine  special  and  general  perturbation 
techniques  to  compensate  for  their  disadvantages  [60]. 

Precise  orbit  prediction  is  of  great  importance  in  the  space  age.  Accurate  positions 
of  operational  satellites  help  us  to  communicate  with  them.  Antennas  are  useless  if  they 
are  not  pointing  to  the  right  direction,  unless  they  are  omnidirectional.  Nonoperational 
satellites  and  debris  pose  threats  to  operational  satellites.  Accurate  orbit  estimation  is  very 
important  for  collision  avoidance  because  wrong  maneuver  decisions  cannot  only  result 
in  the  loss  of  millions  of  dollars  but  also  create  thousands  of  pieces  of  debris.  Moreover, 
accurate  predictions  for  the  reentry  of  returning  space  objects  are  very  important  for  two 
reasons.  First,  space  objects  returning  from  space  might  pose  threats  to  people.  Second, 
space  assets  returning  from  space  might  contain  hardware  or  software  important  to  the 
national  security.  Therefore,  U.S.  Space  Surveillance  Network  (SSN)  detects,  tracks,  and 
catalogs  debris,  operational  and  nonoperational  satellites  every  day.  SSN  has  more  than 
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23,000  objects  orbiting  Earth  in  its  database.  [54]. 

North  American  Aerospace  Defense  Command  (NORAD)  publishes  Two-Line  Ele¬ 
ment  (TLE)  sets  of  some  satellites  to  promote  situational  awareness  in  space.  Special 
General  Perturbations  4  (SGP4)  is  the  current  method  of  NORAD  to  propagate  orbits  from 
raw  observational  data  coming  from  radars.  However,  SGP4  is  accurate  in  the  vicinity  of 
epoch  time.  SGP4  is  based  on  general  perturbation  techniques.  Therefore,  it  suffers  from 
errors  that  emerge  from  omission  of  higher  order  terms  in  the  equations  of  motions,  and  its 
predictions  are  accurate  only  in  kilometers.  SGP4  requires  observational  data  to  make  cor¬ 
rections  once  every  few  days.  However,  it  is  very  fast  in  terms  of  computation  time.  Special 
perturbation  techniques  are  used  if  more  accurate  predictions  are  needed.  For  example,  the 
success  of  a  collision  avoidance  intrinsically  depends  on  accuracy  in  the  positions  of  ob¬ 
jects  in  danger  of  collision.  The  downsides  of  numerical  integration  methods  are  truncation 
errors  in  the  long  run  due  to  the  limited  word  length  of  computers,  high  computational  cost, 
and  difficulty  in  determining  whether  the  resultant  orbit  is  the  one  that  is  desired.  More¬ 
over,  the  breakup  of  satellite  1961-Omicron  in  1961  proved  that  heavily  dependence  on 
numerical  methods  that  requires  high  computational  power  for  the  sake  of  accuracy  is  not 
reliable.  In  1961,  the  breakup  tripled  the  number  of  debris  to  be  tracked.  The  Naval  Ordi¬ 
nance  Research  Calculator  (NORC),  which  was  at  the  Naval  Weapon’s  Laboratory  (NWL) 
in  Dahlgren,  was  insufficient  to  cope  with  the  surveillance  data  flow  and  the  orbital  up¬ 
dates.  Recently  acquired  IBM  7090,  which  was  three  times  faster  than  NORC  in  terms  of 
computational  power,  and  incredible  human  effort  processed  the  data  in  a  few  days.  That 
catastrophic  event  has  been  second  to  none  since  then.  However,  lessons  were  well  learned. 
It  is  crucial  to  have  accurate  orbit  determination  models  that  don’t  require  as  much  compu¬ 
tational  power  as  special  perturbation  techniques  [26].  Therefore,  a  new  method  for  orbit 
prediction,  which  is  as  accurate  as  numerical  methods  and  as  fast  as  analytical  methods,  in 
terms  of  computational  time,  would  be  very  desirable. 
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1.2  Problem  Statement 


This  work  will  answer  the  question  of  whether  best  low  eccentricity  Kolmogorov  Arnold 
Moser  (KAM)  tori  fitted  to  the  orbits,  using  TLE  histories  of  satellites,  produce  more  accu¬ 
rate  predictions  for  low  eccentricity  orbits  than  SGP4  does.  Moreover,  this  effort  presents 
an  approximate  accuracy  analysis  for  the  low  KAM  torus  orbit  prediction  because  both 
SGP4  and  TLE  suffer  from  intrinsic  inaccuracies.  The  actual  accuracy  that  can  be  achieved 
by  this  theory  can  be  analyzed  using  raw  observational  data. 

A  previous  study  by  Frey  showed  that  it  is  possible  to  extract  KAM  torus  basis  fre¬ 
quencies  from  SGP4  and  TLE  sets.  Frey  concluded  that  low  eccentricities  of  Hubble  Space 
Telescope  (HST)  and  Thor  Rocket  Body  caused  difficulties  in  the  construction  process  of 
KAM  torus,  and  the  difficulties  in  the  TLE  curve  fitting  for  the  two  Delta  Rocket  Bod¬ 
ies  are  due  to  air  drag.  Frey  used  a  modified  Laskar  frequency  algorithm,  developed  by 
Wiesel,  to  determine  KAM  torus  basis  frequencies  [28].  Laskar  fequency  analysis  is  used 
to  determine  the  stability  of  orbits  in  dynamical  systems  where  the  energy  is  conserved. 

Because  the  modified  Laskar  frequency  algorithm  cannot  model  low  eccentricity  orbits 
to  the  desired  degree  of  accuracy,  Wiesel  developed  a  new  theory  for  low  eccentricity  Earth 
satellite  motion  to  construct  KAM  torus.  This  work  will  make  use  of  this  new  theory  to 
build  KAM  torus  for  low  eccentricity  orbits.  This  KAM  torus  construction  method  includes 
geopotential,  order  and  degree  20,  and  air  drag  perturbations  [69]. 

1.3  Approach 

This  effort  uses  three  different  main  programs  to  prove  KAM  torus  orbit  prediction 
from  TLE  sets  provides  more  accurate  predictions  in  the  long  run  than  SGP4  and  TLE  sets. 
The  first  main  program  is  a  public  domain  SGP4  program  written  in  C++  by  Vallado  [61]. 
It  outputs  predicted  positions  and  velocities  for  an  orbit  period  and  some  orbital  elements 
by  making  use  of  TLE  sets  of  satellites.  TLE  sets  of  satellites  with  low  eccentricities  are 


4 


obtained  from  www.space-track.org  once  every  3  to  4  days  for  a  period  of  2  months.  NO¬ 
RAD  imports  corrections  to  SGP4  predictions  with  observational  data  coming  from  radars 
once  every  3  to  4  days.  The  second  main  program,  developed  by  Wiesel,  builds  a  low  ec¬ 
centricity  KAM  torus  including  zonal,  sectoral,  and  tesseral  geopotential  perturbations,  and 
air  drag  perturbations.  The  third  main  program  is  a  least  squares  fitting  program  that  fits 
a  low  eccentricity  KAM  torus  to  SGP4  and  TLE  data  by  iteratively  correcting  modal  vari¬ 
ables,  which  are  KAM  torus  coordinates  and  their  linearizations,  and  basis  frequencies  of 
the  KAM  torus,  and  subsequently  outputs  residual  mean  square  values.  The  methodology 
is  based  on  the  concept,  which  is  very  similar  to  that  of  the  least  squares  method,  that  it  is 
possible  to  obtain  more  accurate  results  with  more  data,  and  SGP4  is  more  accurate  when 
nearer  to  the  epoch  time.  The  author  has  developed  some  scripts  to  analyze  all  low  eccen¬ 
tricity  satellites,  which  have  eccentricities  smaller  than  10  2,  and  inclinations  not  close  to 
critical  and  polar  inclination.  There  are  7,938  satellites  and  pieces  of  debris  correspond 
with  the  above  description,  which  can  be  publicly  obtained  from  www.space-track.org. 
However,  the  author  has  pseudo-randomly  chosen  1500  of  these  7,398  near-Earth  objects 
as  his  test  cases.  This  effort  serves  as  a  prelude  to  converting  the  TLE  sets  catalog  from 
SGP4  to  low  eccentricity  KAM  torus  theory. 

1.4  Results 

This  effort  proves  that  the  low  eccentricity  KAM  theory  is  a  better  substitute  for  SGP4. 
The  new  theory  is  five  times  more  accurate  for  the  best  fits,  and  three  times  better  for  mean 
fits.  Most  of  the  low  eccentricity  orbits  of  the  near-Earth  objects  can  be  modeled  by  the  low 
eccentricity  KAM  torus  method.  However,  resonance,  higher  eccentricity,  polar  and  critical 
inclinations,  and  air  drag  issues  should  be  addressed.  The  results  are  certainly  promising. 
Specifically,  the  new  method  can  be  used  for  collision  avoidance  calculations,  and  forma¬ 
tion  flight  applications.  The  new  method  provides  with  a  set  of  numerical  algorithms  that 


5 


may  rival  the  numerical  methods  in  accuracy  and  the  analytic  methods  in  computational 
speed.  This  work  also  presents  a  rigorous  performance  analysis  of  the  theory,  which  will 
be  hopefully  used  as  a  reference  for  future  studies. 

1.5  Overview 

This  document  is  organized  in  five  chapters.  Chapter  II  will  present  previous  studies 
conducted  by  Wiesel  and  his  masters  and  PhD  students  on  KAM  torus  orbit  determination. 
It  will  also  give  an  overview  of  the  history  of  analytical  orbit  modeling  in  the  United  States 
space  surveillance  system,  satellite  dynamics  formulated  in  Hamiltonian  mechanics,  KAM 
theory,  the  geopotential  and  air  drag  as  perturbing  forces,  and  least  squares  fitting.  Chapter 
III  will  discuss  the  methodology  and  the  data  to  be  analyzed.  It  will  also  outline  predicting 
positions  and  velocities  from  TLE  and  SGP4,  building  low  eccentricity  KAM  tori,  and 
fitting  the  low  eccentricity  KAM  tori  to  the  orbits  obtained  from  SGP4  and  TLE.  Chapter  IV 
will  present  the  results  and  limitations  of  KAM  torus  orbit  prediction  method.  It  will  also 
compare  the  accuracy  of  the  new  method  to  that  of  SGP4  and  TLE.  A  general  summary  of 
the  work,  conclusions,  and  recommendations  for  future  studies  will  be  presented  in  Chapter 
V. 
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II.  Background 


This  chapter  provides  background  information  on  current  orbit  determination  methods, 
near-Earth  satellite  motion,  and  KAM  theory.  First,  a  brief  history  of  analytical  orbit  mod¬ 
eling  in  the  United  States  Space  Surveillance  System  is  presented.  Then,  the  importance  of 
the  space  situational  awareness  is  presented.  Next,  Hamiltonian  mechanics,  transformation 
theory,  orbit  perturbations,  and  KAM  theory  are  described.  Then,  the  application  of  KAM 
theory  to  celestial  mechanics  is  presented.  After  that,  the  previous  studies  of  Wiesel  and 
his  masters  and  PhD  students  on  the  applications  of  KAM  theory  to  Earth  satellites  are 
provided  in  chronological  order.  Finally,  a  summary  of  this  chapter  will  be  given. 

2.1  A  Brief  History  of  Analytical  Orbit  Modeling 

After  the  launch  of  Sputnik  in  1957,  United  States  started  to  track  space  objects.  Today, 
SSN  database  has  more  than  23,000  space  objects  orbiting  earth  [54].  Earth  orbiting  objects 
are  cataloged  by  analytical  orbit  models.  The  analytical  methods  have  been  mostly  devel¬ 
oped  by  the  Air  Force  Space  Command  and  the  Naval  Space  Command.  The  development 
and  improvement  of  orbital  models  and  algorithms  span  from  1957  to  today  [26]. 

The  necessity  of  knowing  orbital  positions  of  space  objects  emerged  due  to  military 
concerns  in  1957.  The  Air  Force  used  it  not  to  confuse  a  missile  with  an  object  orbiting 
Earth,  and  the  Navy  used  it  to  warn  the  fleets  against  space  reconnaissance.  The  potential 
benefits  from  artificial  satellites  led  to  great  interest  in  the  accurate  positions  of  the  satellites 
for  the  civilian  world.  Therefore,  the  catalog  for  Earth  orbiting  objects  was  created,  and  it 
is  still  in  use  today  [26] . 

The  first  formal  catalog  was  created  at  the  National  Space  Surveillance  Control  Center 
(NSSCC).  The  control  center  was  located  in  Bedford,  Massachusetts.  The  observational 
data  were  provided  by  150  different  sites.  There  were  4  different  types  of  observations, 
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which  were  categorized  by  their  source  and  content.  Table  1  shows  details  of  different 
observational  data  types  used  from  the  year  1957  to  1963  [26]. 

Table  1.  Observational  Data  Types  [26] 


Observation  Type 

Content 

Source 

Type  1 

2  angles  and  slant  range 

Radars 

Type  2 

2  angles 

Baker-Nunn  cameras, 
telescopes,  binoculars, 
visual  sightings 

Type  3 

Azimuth 

Direction  finders 

Type  4 

Time  of  closest  approach  (Doppler) 

Radars,  radio  receivers 
(for  transmitting  satellites) 

Observations  were  processed  by  an  IBM-709  computer  to  update  orbital  elements  in 
the  catalog,  which  were  some  type  of  mean  orbital  elements.  Then,  updated  orbital  data 
were  processed  by  another  program  to  yield  three  different  products,  which  were  used  to 
make  observations  for  the  next  time  by  ground  sites.  The  first  product,  called  bulletin,  was 
the  precursor  of  the  TLE  sets.  The  most  important  components  of  a  bulletin  were  longi¬ 
tude,  time  of  ascending  pass,  and  revolution  number  of  a  satellite.  These  were  generated 
approximately  a  week  in  advance.  The  second  product  was  used  to  calculate  all  satellite 
passes,  which  could  be  seen  from  an  observing  site.  The  third  product  was  very  similar 
to  the  second  product.  However,  the  third  one  was  used  by  the  U.S.  Navy  and  U.S.  Army 
observation  fences,  which  were  composed  of  a  number  of  radars.  Therefore,  the  third  prod¬ 
uct  was  used  to  yield  the  intersection  point  of  orbital  plane  of  a  satellite  and  vertical  plane 
spanned  by  radar  beam  instead  of  some  look  angles  and  ranges  [26]. 

In  1961,  U.S.  Navy  constructed  the  Naval  Space  Surveillance  System  (NAVSPASUR), 
which  detected  and  cataloged  Earth  orbiting  objects  mostly  without  human  intervention. 
The  processing  unit  for  NAVSPASUR  was  located  at  Dahlgren,  Virginia.  The  system  was 
a  continuous  wave  multi-static  interferometer,  which  observed  LEO  satellites  4  to  6  times  a 
day.  There  were  3  transmitters  and  6  receivers,  which  spanned  from  San  Diego,  California 
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to  Savannah,  Georgia  as  9  individual  sites  located  in  the  Southern  States.  The  system 
was  called  “Fence”.  The  command  of  the  Fence,  which  was  responsible  for  40%  of  all 
observations  made  by  the  Air  Force,  was  passed  to  U.S.  Air  Force  in  2004.  The  U.S.  Air 
Force  shut  down  the  system  in  October  2013.  A  new  fence  is  planned  to  be  operational  in 
2017  in  the  Marshall  Islands  [19]. 

Special  perturbation  techniques  were  thought  to  be  used  in  the  Fence  until  it  was  dis¬ 
covered  that  it  took  the  NORC  at  the  NWL  several  hours  to  update  one  orbit.  In  1961, 
the  catalog  was  transferred  to  an  IBM-7090  computer,  which  reduced  the  time  required 
to  update  one  orbit  to  approximately  1  minute.  That  was  very  promising  for  the  imple¬ 
mentation  of  highly  numerical  methods,  which  yield  more  accurate  predictions.  However, 
the  breakup  of  1961-Omicron  tripled  the  number  of  objects  to  be  tracked.  The  number 
of  detectable  objects  exceeded  the  capability  of  the  NORC.  The  IBM-7090  and  incredible 
human  effort  processed  the  data  in  a  few  days.  Although  special  perturbation  techniques 
provide  with  more  accurate  predictions,  the  beakup  of  satellite  1961-Omicron  proved  that 
accurate  analytical  models  are  needed  in  case  computers  fail  to  process  the  data  [26]. 

In  1959,  Dirk  Brouwer  and  Yoshide  Kozai  provided  two  different  solutions  to  the  same 
problem,  which  was  Earth  satellite  motion  under  the  influence  of  the  zonal  harmonics  J2, 
J3,  J4,  and  J5  [9,  38].  Modern  orbit  prediction  methods  used  in  the  U.S.  Space  Surveillance 
System  include  fundementals  of  the  solutions  provided  by  Brouwer  and  Kozai.  In  1961, 
Brouwer  and  Gen-Ichiro  Hori  added  atmospheric  drag  effect  to  the  1959  Brouwer  solution 
[10].  However,  the  atmospheric  model  was  computationally  heavy  for  the  computers  of 
the  time  due  to  the  series  expansions  in  the  scale  height.  In  1960,  a  group  of  people  under 
the  project  SPACETRACK  developed  an  atmospheric  density  modeling,  which  solved  the 
series  expansion  problem  for  artificial  satellites.  Max  Lane  developed  an  analytic  orbit 
model  based  on  the  atmospheric  density  modeling  of  the  SPACETRACK  in  1965  [40]. 
Lane  and  Cranford  improved  the  analytic  orbit  model  by  implementing  analytic  density 
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model  instead  of  empirical  one  in  1969  [42].  In  1963,  Lyddane  contributed  to  the  analytic 
orbit  modeling  by  solving  small  divisors  of  the  eccentricity  and  sine  of  the  inclination 
by  defining  the  perturbation  theory  in  terms  of  Poincare  variables  instead  of  Delaunay 
variables  [49]. 

Conversion  from  theory  to  practice  for  analytical  orbit  models  took  place  in  Dahlgren, 
Virginia,  and  Colorado  Springs,  Colorado.  NAVSPASUR  implemented  the  1959  Brouwer 
solution  with  Lyddane’s  contribution,  which  was  known  as  PPT  (Position  and  Partials  as  a 
function  of  Time),  on  an  IBM-7090  computer  in  1964.  PPT  adopted  semi-empirical  drag 
model  developed  by  Richard  H.  Smith  because  the  atmospheric  model  of  Brouwer  and  Hori 
was  computationally  heavy  for  the  computers  of  the  time,  and  it  was  too  early  for  the  Lane’s 
analytic  orbit  model.  The  change  in  the  eccentricity  could  be  solved  using  Equation  1  and 
Equation  2. 


eD  =  e0(l-e02)  — 

do 


(1) 


4  a0  d() 

a”  =  _3^(T) 


(2) 


Numerically,  mean  motion  term  in  PPT  is  similar  to  Kozai’s  mean  motion,  which  has  the 
zonal  secular  perturbation  rate  of  mean  anomaly  drived  by  Brouwer. 

NSSCC  was  moved  to  Colorado  Springs,  and  named  as  Space  Detection  and  Tracking 
System  (SPADATS).  SPADATS  is  the  other  location  for  the  implementation  of  the  ana¬ 
lytic  orbital  theory,  which  was  called  SGP  (Simplified  General  Perturbations),  developed 
by  Brouwer  and  Kozai.  Arsenault,  Chaffee,  and  Kuhlman  avoided  small  divisors  of  the 
eccentricity  and  sine  of  the  inclination  in  SGP  by  defining  the  solution  in  non-singular 
terms  and  keeping  only  the  most  important  ones  in  the  model.  They  included  long  and 
short-period  terms  without  the  eccentricity  effect  from  the  Brouwer  solution  and  the  con- 
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vention  that  relates  mean  motion  to  semi-major  axis  from  the  Kozai  solution  [33].  The  air 
drag  effect  used  in  SGP  was  very  similar  to  Smith’s  semi-empirical  model.  In  1969,  Lane 
and  Cranford  developed  an  analytic  atmospheric  model  for  artificial  satellites  [42].  How¬ 
ever,  computers  of  the  time  were  inadequate  to  cope  with  calculations  required  for  the  new 
model.  Therefore,  the  most  important  terms  of  the  atmospheric  model  were  implemented 
to  the  analytical  orbit  model,  which  is  now  called  SGP4.  Today,  some  flavor  of  SGP4  is 
used  by  NORAD  [26]. 

Solar  and  Lunar  gravitation,  and  resonance  in  the  tesseral  harmonics  of  Earth  became 
important  when  the  first  satellite  with  a  period  of  12  hours  was  launched  in  1965.  Bruce 
Bowman  developed  a  model  which  includes  third  body  effects,  and  resonance  in  the  tesseral 
harmonics  of  Earth  in  1967  [7].  In  1977,  Dick  Hujsak  developed  another  model  which  in¬ 
cludes  everything  from  Bowman’s  work.  Moreover,  his  first-order  model  could  be  applied 
to  geosynchronous  satellites  [32].  Hujsak’s  model  was  combined  with  the  SGP4  model, 
which  is  still  in  use  today.  In  1997,  third  body  effects,  and  resonance  in  the  tesseral  har¬ 
monics  of  Earth  were  adopted  from  SGP4  by  the  Naval  Space  Command.  The  new  method, 
which  is  still  in  use  today,  is  called  PPT3  [26]. 

2.2  Space  Situational  Awareness 

In  less  than  six  years  from  the  launch  of  Sputnik  1  to  1963,  616  man-made  objects 
accumulated  in  low  Earth  orbit  (LEO).  Those  man-made  objects  were  76  payloads,  35 
rocket  bodies,  and  227  mission  related  debris.  The  U.S.  Ablestar  upper  stage  explosion 
created  184  pieces  of  objects  which  were  half  of  all  the  objects  cataloged  by  SSN  in  1963, 
and  60%  of  all  debris  from  this  explosion  is  still  orbiting  Earth  [54].  Today,  23,000  objects 
are  cataloged  by  SSN  and  the  majority  of  these  objects  are  not  operational  satellites  but 
debris  and  nonoperational  satellites,  which  is  approximately  95%  of  total  number  of  objects 
being  cataloged  [59]. 
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Figure  1.  Cataloged  objects  by  SSN  in  1963  (left)  and  in  2013  (right)  [54] 

The  rate  of  growth  in  the  number  of  debri  objects  in  less  than  six  years  was  noticed 
by  Ernest  Peterkin,  Head  of  the  Systems  Section  of  the  Operational  Research  Branch  at 
the  U.S.  Naval  Research  Laboratory.  He  was  one  of  the  first  to  calculate  the  number  of 
objects  per  unit  volume  in  LEO  in  Lebruary  1963.  He  derived  population  growth  functions 
for  the  objects  in  LEO  using  estimates  for  launch  rates,  lifetimes  of  satellites,  and  satellite 
collisions.  Peterkin’s  conclusions  showed  that  there  would  be  16,500  objects  in  LEO  at 
the  beginning  of  2013.  His  predictions  are  very  close  to  the  number  of  objects  cataloged 
by  SSN  on  1  January  2013.  Peterkin  concluded  that  space  surveillance  systems  would 
have  difficulty  coping  with  the  rapid  growth  in  the  number  of  objects  orbiting  Earth  in 
the  future.  He  stated  that  it  is  important  to  improve  detection  and  data  processing  methods 
used  in  space  surveillance  systems  to  manage  the  increase  in  the  number  of  objects  orbiting 
Earth  [54].  This  effort  is  trying  to  find  an  accurate  and  a  fast  orbit  prediction  method 
which  can  help  space  surveillance  systems  to  cope  with  the  proliferation  of  objects  orbiting 
Earth.  Low  eccentricity  KAM  theory  yields  more  accurate  orbit  predictions  than  SGP4.  In 
addition,  this  work  shows  that  it  is  possible  to  extract  valuable  data,  which  are  used  by  the 
low  eccentricity  KAM  torus  orbit  prediction  method  to  make  corrections  to  low  eccentricity 
KAM  tori,  from  inaccurate  TLEs  in  the  long  run. 

The  Joint  Space  Operations  Center  (JSpOC),  located  at  Vandenberg,  is  responsible  for 
detecting,  tracking,  and  cataloging  all  man-made  objects  orbiting  Earth.  JSpOC  uses  SSN 
which  has  30  space  surveillance  radars  and  optical  telescopes  all  over  the  world.  The  radar 
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Figure  2.  Worldwide  network  of  space  surveillance  sensors  of  SSN  [55] 

sensors  are  dedicated  to  make  observations  for  near-Earth  objects,  which  are  below  6000 
kilometers  altitude.  The  electro-optic  sensors  are  dedicated  to  make  observations  for  deep- 
space  objects,  which  are  above  6000  kilometers  altitude.  Azimuth,  elevation,  range,  and 
range  rate  are  typical  observation  types  obtained  from  a  radar  site.  Azimuth  and  elevation 
or  right  ascension  and  declination  are  observations  provided  by  an  optical  sensor.  Observa¬ 
tions  made  by  SSN  include  time  tags  [37].  Figure  2  shows  locations  of  worldwide  network 
of  30  space  surveillance  sensors.  Observational  data  obtained  by  these  sensors  are  fitted  to 
the  trajectories  of  orbiting  objects,  and  the  satellite  catalog  is  continuously  updated.  SSN 
makes  380,000  to  420,000  observations  every  day.  SSN  implements  a  predictive  surveil¬ 
lance  method  because  all  objects  orbiting  Earth  cannot  be  covered  by  existing  sensors, 
which  requires  efficient  analytic  orbit  models.  These  observations  are  important  for  opera¬ 
tional  satellites  because  collision  analyses  are  made  using  orbit  predictions  in  the  satellite 
catalog  [59]. 

This  effort  is  beneficial  for  satellite  communities  which  use  NORAD  TLE  for  propa¬ 
gating  orbits  of  the  objects  that  pose  threat  to  their  space  assets  because  low  eccentricity 
KAM  theory  can  yield  position  and  velocity  data  to  create  TLEs  before  NORAD  publishes 
them.  NORAD  publishes  a  new  element  set  when  the  difference  of  positions  predicted  by 
the  current  and  the  new  element  sets  exceed  5  km  with  a  90%  confidence  interval.  LEO 
satellites  require  more  frequent  updates  due  to  the  atmospheric  drag.  Maneuverable  satel- 
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lites  also  need  more  frequent  updates  because  orbital  maneuvers  cannot  be  modeled  by 
SGP4  [37].  This  work  uses  Vallado’s  public  domain  SGP4  with  differential  corrections  to 
determine  the  accuracy  of  the  orbit  predictions  yielded  by  low  eccentricity  KAM  theory 
because  JSpOC  version  of  SGP4  has  no  technical  details  and  codes  available  to  public. 

2.3  Orbit  Prediction  using  SGP4  and  TLE 

TLEs  are  some  type  of  mean  elements  which  are  averaged,  and  smoothly  changing  over 
time  or  angle.  Some  period  of  time,  mean  anomaly,  true  anomaly,  and  eccentric  anomaly 
can  be  used  to  average  the  elements  [60].  TLEs  are  used  by  NORAD,  and  the  majority  of 
the  satellite  community.  NORAD  TLE  is  the  default  input  for  most  of  the  commercialized 
satellite  ground  antenna  control  system  and  orbit  analysis  programs  [46].  TLE  is  gener¬ 
ated  from  observational  data  flow  from  SSN  using  one  of  five  orbit  prediction  methods 
of  NORAD.  The  first  one  is  SGP  which  was  developed  by  Hilton  and  Kuhlman  [31],  and 
simplified  by  the  work  of  Kozai  in  terms  of  gravitational  and  air  drag  effects  [38].  SGP 
is  used  for  satellites  with  a  period  less  than  225  minutes.  The  second  one  is  SGP4  which 
was  developed  by  Cranford  in  1970  [41],  and  simplified  by  adopting  analytic  atmospheric 
model  of  Lane  and  Cranford  instead  of  empirical  atmospheric  model  in  1969  [42].  Table  2 
shows  parameters  required  to  initialize  SGP4.  All  of  the  parameters  but  mean  motion  are 
mean  orbital  elements  defined  by  Brouwer  [9].  Brouwer’s  mean  elements  include  only  the 
effects  of  the  zonal  harmonics  which  are  J2,  J3,  J4,  and  J5.  SGP4  includes  orbital  mean 
elements  and  their  linearizations  defined  in  series  expansions.  Many  assumptions  are  made 
for  SGP4.  For  example,  Earth’s  gravitational  potential  includes  only  a  few  zonal  terms, 
the  atmospheric  model  is  a  static  model  with  an  assumption  of  exponential  decay,  and  the 
third-body  mass  and  resonance  effects  are  partially  added  to  the  model  [62].  Mean  motion 
in  the  Table  2  includes  only  short-period  oscillations  similar  to  Kozai’s  mean  motion  [38]. 
Long-periodic  oscillations  are  masked  by  air  drag,  which  gradually  increases  mean  motion 
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Figure  3.  TLE  Example  for  NOAA-6  Weather  Satellite 


over  time,  for  LEO  satellites  [28]. 

Table  2.  Parameters  to  Initialize  SGP4  Propagation 


Symbols 

Description 

to 

Epoch  time 

na 

Mean  motion  at  epoch 

ea 

Eccentricity  at  epoch 

io 

Inclination  at  epoch 

(Oo 

Argument  of  perigee  at  epoch 

o0 

Right  ascension  of  the  ascending  node  at  epoch 

M0 

Mean  anomaly  at  epoch 

B* 

Atmospheric  drag  coefficient 

The  third  one  is  SDP4  which  is  developed  for  orbits  with  a  period  more  than  225  min¬ 
utes  by  Hujsak  in  1979  [32].  SDP4  includes  third-body,  sectoral,  and  tesseral  effects, 
whereas  SGP4  doesn’t.  The  fourth  one  is  SGP8  which  is  applicable  to  near-Earth  satellites. 
SGP8  depends  on  the  model  of  Lane  and  Cranford  for  gravitational  and  air  drag  effects 
[42].  The  integration  of  differential  equations  are  treated  in  a  different  way  for  SGP8  than 
SGP4.  The  fifth  one  is  SDP8  which  is  a  deep-space  model.  SGP8/SDP8  were  introduced  to 
decrease  the  effects  of  the  deficiencies  of  SGP4/SDP4  associated  with  reentry  and  orbital 
decay  [62].  A  TLE  should  be  used  with  one  of  these  models  because  periodic  variations 
are  removed  in  a  way  that  these  five  models  can  compensate  for  them.  Figure  3  shows  an 
example  of  a  NORAD  TLE  with  descriptions  of  each  elements. 

Orbit  prediction  methods  are  divided  into  three  different  catagories:  special  perturba- 
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tion  techniques,  general  perturbation  techniques  and  semianalytical  techniques.  Special 
perturbation  techniques  use  numerical  methods  to  make  predictions.  However,  numerical 
methods  are  prone  to  errors  due  to  limited  word  length  of  computers.  Moreover,  long  time 
steps  for  integrations  result  in  round-off  errors.  Special  perturbation  techniques  cannot  be 
generalized  to  other  orbits  because  the  predictions  are  special  to  the  orbit  that  is  integrated. 
Special  perturbation  techniques  provide  predictions  accurate  in  meters.  However,  they  re¬ 
quire  high  computational  power.  All  perturbations  are  calculated  for  every  point  in  time 
in  special  perturbation  techniques.  General  perturbation  techniques  use  closed  form  solu¬ 
tions  to  equations  of  motions  of  orbits.  The  complexity  of  equations  of  motions  requires 
omission  of  higher  order  terms  in  the  equations.  The  accuracy  of  the  predictions  degrades 
rapidly  due  to  these  omissions.  In  addition,  orbit  predictions  made  by  general  perturbation 
techniques  are  not  as  accurate  as  those  provided  by  special  perturbation  techniques.  How¬ 
ever,  general  perturbation  techniques  provide  orbit  predictions  without  calculating  every 
intermediate  points  to  find  the  desired  point  in  the  future  and  take  an  average  of  all  pertur¬ 
bations  over  time  as  one  parameter  which  can  be  integrated  into  the  equations.  Therefore, 
they  are  very  fast  in  terms  of  computation  time.  SGP4  orbit  determination  method  is  such  a 
general  perturbation  technique.  Semianalytical  techniques  combine  advantages  of  analyti¬ 
cal  methods  and  numerical  methods. 

Equation  3  shows  the  relationship  between  the  SGP4  model  and  the  state  vector  of  a 
satellite. 

y(t)=fsGP4(x0,B*,t),  (3) 

where  /sgpa  is  the  SGP4  model,  y(t)  is  the  velocity  and  position  vectors  at  time  t,  x0  is 
the  mean  orbital  elements,  which  are  shown  in  the  Table  2  excluding  ta  and  B *,  at  the 
epoch  time.  The  epoch  time  is  represented  as  t0  [46].  Equation  4  shows  how  the  ballistic 
coefficient,  B ,  and  the  atmospheric  reference  reference  density  p0  are  related  to  NORAD 
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SGP4  air  drag  term  B*  [52]. 


B*  =  l-BPo  (4) 

Equation  5  can  also  be  used  to  produce  the  NORAD  TLE  given  osculating  elements  and 
a  good  approximation  for  the  NORAD  TLE.  Equation  5  shows  the  relationship  between  os¬ 
culating  elements  and  the  NORAD  mean  elements.  B*  is  missing  in  the  Equation  5  because 
SGP4  can’t  model  the  air  drag  term.  Therefore,  B*  needs  to  be  calculated  separately  [46]. 
B*  is  important  for  the  accuracy  of  the  NORAD  TLE.  The  omission  of  B*  value  of  10  4 
could  worsen  the  accuracy  of  SGP4  twice  as  much  within  approximately  3  days  [45]. 

yi=  fsGP4i(Xi,...,X6),  (5) 

where  i=(  1  -  ,6),  and  y,-  are  osculating  elements,  jq  are  the  NORAD  mean  elements.  Equa¬ 

tion  5  can  be  expanded  around  xf-,  which  are  good  approximations  for  the  NORAD  mean 
elements,  using  Taylor  series  expansion.  Equation  6  and  Equation  7  show  the  linear  system 
of  equations  after  taylor  series  expansion  [45]. 

Ayi  =  MAxi ,  (6) 

dfsGP4x  dfscP4l  dfsGP4l 

dtf  dx°2  dx°6 

m=  \  (7) 

dfsGP46  dfsGP4b  dfsGP46 

drf  dx%  dx^ 

where  Av,-  =  y;  —y“,  and  At,  =  x,  —  x“ .  Partial  derivatives  in  matrix  M  are  hard  to  solve  due 
to  the  coupled  nature  of  the  SGP4  model.  Newton’s  quotient  can  be  used  to  find  the  partial 
derivatives  [25].  Equation  8  shows  one  of  the  partial  derivatives  in  the  difference  quotient 
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notation  [45]. 


dfsGPA,  =  fsGPlirf  +  M^2,  •  •  •  >-*e)  -  fsGPA (A  i ; %2 > ' ' '  (8) 

c)xax  Ax“  ’  1  J 

where  Ar"  can  be  some  small  increment  of  percentage  of  x“ .  The  NORAD  mean  elements 
can  be  approximated  using  a  process  called  differential  corrections.  The  iterative  process 
is  stopped  when  the  difference  between  the  given  osculating  elements  and  derived  ones, 
which  are  (v,  —  y°),  are  small  enough  for  convergence.  Equation  8  shows  the  iterative 
process  to  approximate  the  NORAD  mean  elements  [45]. 

*?+1=*?  +  EV(yf-:v?),  (9) 

7=1 

where  i=(l,  -  •  •  ,6). 

There  are  inevitable  deficiencies  in  the  SGP4  theory.  Because  SGP4  is  an  analytic  or¬ 
bit  model,  it  ignores  higher  order  terms  in  the  equations  of  motion.  Moreover,  sectoral 
and  tesseral  gravity  field  perturbations  are  not  included  in  the  model.  The  SGP4  model 
includes  only  the  effects  of  the  zonal  harmonics  which  are  J2,  J3,  J4,and  J5.  Therefore,  the 
SGP4  model  yields  position  errors  of  2  km  at  epoch  time  [52].  The  velocity  predictions 
made  by  the  NORAD  TLE  and  the  SGP4  are  less  accurate  than  the  position  predictions  be¬ 
cause  the  rates  of  change  in  the  orbital  parameters,  which  are  not  accurate  due  to  inevitable 
deficiencies  in  the  SGP4  theory,  are  used  to  calculate  velocities  [28]. 

The  accuracy  of  the  Vallado’s  SGP4  with  differential  corrections  for  near-Earth  objects 
are  on  the  order  of  magnitude  of  100  km.  The  prediction  accuracy  worsens  with  decreasing 
altitude,  and  increasing  eccentricity  and  solar  flux.  Table  3  and  Table  4  show  the  SGP4 
orbit  prediction  accuracy  of  the  near-Earth  objects  in  terms  of  the  root  mean  square  values 
with  respect  to  two  different  solar  flux  values  [24]. 
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Table  3.  The  Orbit  Prediction  Accuracy  of  the  Near-Earth  Objects  (F]0  =  100)[24] 


Altitude(km) 

Test  Cases 

lday(km) 

3days(km) 

7days(km) 

15days(km) 

h<400 

57 

4 

10 

60 

300 

400<h<600 

168 

3 

10 

50 

100 

600<h<1200 

267 

3 

10 

20 

50 

1200<h<7000 

90 

2 

10 

10 

20 

Table  4.  The  Orbit  Prediction  Accuracy  of  the  Near-Earth  Objects  ( F\  0  =  200)[24] 


Altitude(km) 

Test  Cases 

lday(km) 

3days(km) 

7days(km) 

15days(km) 

h<400 

57 

10 

40 

300 

1000 

400<h<600 

168 

7 

30 

200 

400 

600<h<1200 

267 

6 

15 

70 

100 

1200<h<7000 

90 

2 

10 

10 

20 

2.4  Reference  Frames 

Reference  frames  are  important  for  classical  and  analytic  mechanics.  All  equations 
of  motions  are  written  in  some  frame  of  reference.  The  first  step  in  solving  a  dynamical 
problem  is  defining  the  reference  frame.  There  are  two  different  types  of  reference  frames, 
which  are  inertial  and  rotating  frames.  There  is  no  solution  for  any  dynamical  problem 
in  classical  mechanics  unless  an  inertial  frame  is  defined.  This  effort  formulate  satellite 


^  =  e 
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Table  5.  Reference  Frames 


Acronym 

Reference 

Frame 

Type 

Origin 

Ist  axis 

2nd  axis 

3  rd  axis 

RAN 

Radial, 

Along- 

Track,and 

Normal 

Rotating 

frame 

Space 

object 

CoM 

Toward 

radial 

direction 

(U) 

Normal 

to  orbital 

frame 

(W) 

Along 

velocity 

direction 

(V) 

ECI 

Earth- 

Centered- 

Inertial 

frame 

Inertial 

frame 

Earth 

CoM 

Toward 

vernal 

equinox 

(7) 

Completes 
right- 
handed 
frame  (£) 

Earth’s 

axis  of 
rotation 
(0 

ECEF 

Earth- 

Centered- 

Earth- 

Fixed 

frame 

Rotating 

frame 

Earth 

CoM 

Toward 
prime 
meridian 
in  equa¬ 
torial 
plane  (x) 

Completes 
right- 
handed 
frame  (y) 

Earth’s 

axis  of 
rotation 

(z) 

ECNF 

Earth- 

Centered- 

Node- 

Fixed 

frame 

Rotating 

frame 

Earth 

CoM 

Toward 
ascend¬ 
ing  node 
in  equa¬ 
torial 
plane  (£) 

Completes 

right- 

handed 

frame 

(V) 

Earth’s 

axis  of 
rotation 
(0 

motion  using  4  different  reference  frames,  which  are  defined  in  Table  5  and  Figure  4. 

Reference  frames  in  Table  5  are  helpful  in  different  ways  for  this  effort.  RAN  reference 
frame  is  used  to  make  residuals  between  SGP4  and  low  eccentricity  KAM  theory  predic¬ 
tions  meaningful  because  it  is  hard  to  analyze  the  residuals  in  ECI  frame.  ECNF  reference 
frame  is  important  because  orbits  are  not  periodic  in  ECI  frame  although  they  are  in  ECNF 
frame.  Periodic  orbits  are  not  periodic  in  ECI  frame  because  orbital  plane  regresses  due  to 
the  zonal  potential.  Moreover,  the  Keplerian  frequency  is  defined  in  a  frame  which  is  tied 
to  the  node  of  the  orbit.  Thus,  it  can  simply  be  generalized  to  ECI  and  ECEF  frame  because 
the  zonal  potential  is  symmetric  about  the  earth’s  polar  axis.  Greenwich  referenced  ECEF 
frame  simplifies  calculations  because  the  gravitational  field  potential  is  constant  in  ECEF 
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frame. 


2.5  Analytical  Dynamics 

Equations  of  motions  are  formulated  using  momentum  and  force  in  Newtonian  me¬ 
chanics.  Objects  and  forces  applied  on  them  are  considered  separately.  For  Newtonian 
mechanics,  constraint  forces,  which  are  forces  that  do  no  work,  may  be  needed  to  solve 
dynamical  problems.  However,  analytical  dynamics  considers  a  dynamical  system  as  a 
whole  not  separately.  For  analytical  dynamics,  equation  of  motions  are  formulated  using 
kinetic  energy  and  work,  which  are  both  scalar  quantities.  Therefore,  constraint  forces  are 
not  needed  to  solve  dynamical  problems.  Dynamical  systems  in  analytical  dynamics  are 
defined  in  generalized  coordinates,  which  are  finite,  continuous,  and  differentiable  with 
respect  to  time.  In  addition,  generalized  coordinates  aren’t  restricted  to  be  physical  quan¬ 
tities.  Number  of  generalized  coordinates  for  a  system  equals  to  the  number  of  degrees  of 
freedom  that  the  system  has  [51].  Many  complex  dynamical  systems  can  be  solved  without 
much  effort  using  analytical  dynamics.  Figure  5  shows  the  evolution  of  transition  from 
Newtonian  to  Fagrangian  mechanics. 

Principle  of  Virtual  Work 

y 

D'  Alembert  Principle 

y 

Hamilton's  Law 

y 

Hamilton's  Principle  (Hamilton's  Extended  Principle) 

y 

Lagrange's  Equations 

Figure  5.  Evolution  of  the  Transition  from  Newtonian  to  Lagrangian  Mechanics 

Principle  of  virtual  work  is  based  on  static  equilibrium  of  a  system.  It  is  the  first  varia¬ 
tional  principle  in  mechanics  [51].  Constraint  forces  are  perpendicular  to  virtual  displace- 
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ments.  However,  formulation  of  equations  of  motions  relies  on  forces.  Equation  10  and 
Equation  11  show  the  virtual  work  of  the  entire  system  equals  to  zero. 

_  N  N 

5  W  =  £  f)  •  8n  +  £  F-  ■  8n  =  0,  (10) 

i=l  i=l 

N 

£  F'i  ■  8rt  =  0,  (11) 

i=  1 

where  F,  are  external  forces,  and  F[  are  constraint  forces  applied  on  the  system,  dr,  are 
virtual  displacements. 

D’  Alembert’s  principle  is  an  extension  of  principle  of  virtual  work  to  dynamics.  D’ 
Alembert’s  principle  treats  dynamical  problems  as  if  they  were  statical.  Equation  12  shows 
the  most  general  formulation  of  dynamics  [51]. 

N 

£  (Pi  ~  Pi)  8ri  =  0,  (12) 

1=1 

where  pi  is  rate  of  change  of  the  momentum,  which  is  also  called  inertia  force. 

Hamilton’s  law  is  important  because  it  depends  on  kinetic  and  potential  energy  instead 
of  force  and  momentum.  Both  principal  of  virtual  work  and  D’  Alembert’s  principle  are 
similar  to  Newtonian  mechanics  in  formulating  equations  of  motion  because  they  all  de¬ 
pend  on  force  and  momentum.  The  first  step  in  transition  from  D’  Alembert’s  principle 
to  Hamilton’s  law  is  integrating  D’  Alembert’s  principal  equation  over  a  finite  period  of 
time,  and  expressing  momentum  in  terms  of  kinetic  energy.  The  derivation  of  Hamilton’s 
law  introduces  Lagrangian,  which  is  the  difference  between  kinetic  and  potential  energy. 
Equation  13  represents  Hamilton’s  law  of  varying  action.  Hamilton’s  law  of  varying  action 
is  valid  for  non-linear,  non-conservative,  non-stationary  dynamic  systems.  Variations  of 
the  generalized  coordinates  Sqj  aren’t  necessarily  zero  because  the  system  state  qj  isn’t 
always  known  at  time  t\  and  h,  which  shows  the  system  is  non-stationary.  Variation  of  any 
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constant  value  equals  to  zero  [39]. 


(13) 


L—T  —  V,  (14) 

where  L  is  Lagrangian,  T  is  kinetic  energy, and  V  is  potential  energy. 

Hamilton’s  principle  is  an  extension  of  Hamilton’s  law.  For  Hamilton’s  principle,  the 
generalized  coordinates  Sqj  are  assumed  to  be  known  at  end  points  t\  and  t2-  Hamilton’s 
principle  is  applicable  when  there  are  no  non-conservative  forces,  whereas  extended  Hamil¬ 
ton’s  principle  is  applicable  when  there  are  non-conservative  forces  in  the  system.  Equa¬ 
tion  15  shows  Hamilton’s  principle,  and  Equation  16  shows  extended  Hamilton’s  principle 
[39].  Hamilton’s  principle  is  a  variational  principle  which  reduces  a  dynamical  problem 
to  a  scalar  integral  independent  of  coordinates,  which  describes  Lagrangian.  Equations  of 
motion  are  obtained  by  figuring  out  the  conditions  which  make  the  scalar  integral  station¬ 
ary.  Figure  6  represents  true  path  and  varied  path,  which  the  true  path  coincide  with  two 
end  points  t\  and  h-  The  varied  path  is  formed  when  the  virtual  displacements  is  applied, 
which  by  definition,  has  no  change  in  5r,  [51]. 


0 


(15) 


ft 2 


8Ldt  +  /  8Wdt  =  0 


(16) 


Lagrange  equations  are  n  number  of  coupled  second-order  differential  equations,  which 
has  n  number  of  generalized  coordinates.  Lagrange  equations  are  derived  from  Hamilton’s 
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Figure  6.  True  Path  and  Varied  Path 


law  of  varying  action  as  [39]: 


N 


dT 


8L(qi,qi,t)dt+  8W-Y(^8qj) 
Jn  dqj 


t2 


=  0 


(17) 


After  calculating  the  variation  of  the  Lagrangian, 


N 


d 2 ,  dL 


dL 


rt 2 . 


N 


dT 


Y  (^8qi+j-Sqi)dt+  SW-J^ij—Sqj) 

Jh  dqi  dq\  Jh  -=1  dqj 


=  0 


(18) 


After  integrating  Lagrangian  terms  by  parts,  and  cancelling  the  trailing  term,  which  shows 
the  system  could  be  non- stationary, 


if 


d  dL  dL 

dt  dqi  dqj  ' 


Sqjdt  =  0, 


(19) 


where  Qj  is  generalized  force.  The  integrand  can  be  equated  to  zero  because  the  system 
can  be  considered  stationary.  Equation  21  shows  Lagrange’s  equations  when  there  is  no 
non-conservative  force  in  the  system  because  virtual  work  equals  to  zero. 


ddL  dL 
dt  dqi  dqi  1 


(20) 
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0 


(21) 


J  ^  r)L  .  dL 
dt  dqi  dqi 


2.6  Hamilton’s  Equations 


The  Hamiltonian  formulation  is  an  alternative  to  the  Lagrangian  formulation  [16].  They 
both  have  the  same  physics.  The  Lagrangian  formulation  has  n  number  of  second-order 
differential  equations,  whereas  the  Hamiltonian  formulation  has  2n  number  of  first-order 
differential  equations.  The  Lagrange’s  equations  describe  mechanics  in  terms  of  general¬ 
ized  coordinates  and  velocities,  whereas  the  Hamilton’s  equations  rely  on  the  generalized 
coordinates,  qt,  and  momenta,  pi  [29].  Equation  22  shows  the  equation  for  the  generalized 
momenta.  The  generalized  momenta  can  be  described  as  linear  functions  of  the  general¬ 
ized  velocities.  Conversely,  the  generalized  velocities  can  be  shown  to  be  linear  functions 
of  the  generalized  momenta  [51].  The  Legendre  transform  of  the  Lagrangian  yields  the 
Hamiltonian.  Equation  23  shows  the  Hamiltonian. 


Pi  = 


dL(qi,qi,t) 

dcp 


(22) 


H(p,q,t)  =  Y,qiPi^L(q,qd),  (23) 

i 

where  H(p,qp)  is  Hamiltonian,  which  can  fully  describe  the  motion  [51]. 

The  Hamilton’s  formalism  are  more  powerful  than  Lagrange  formalism.  The  apparent 
advantage  of  the  Hamilton’s  equations  over  Lagrange  equations  is  that  time  derivatives  of 
the  variables  are  on  the  left  side  of  the  equations,  which  may  help  to  figure  out  first  integrals 
of  motion.  Moreover,  the  Hamilton’s  equations  are  favorable  in  representing  the  solution 
in  geometrically  because  2n-dimensional  representation  is  advantageous  in  representing 
not  only  one  path  but  totality  of  paths,  which  are  all  possible  solutions.  In  Hamiltonian 
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phase  space,  the  velocity  of  every  point  can  be  determined  uniquely,  whereas  the  motion 
can  start  in  any  direction,  and  paths  may  intersect  in  the  Lagrangian  configuration  space, 
which  makes  totality  of  paths  confusing  [51].  The  canonical  equations  of  Hamiltonian  can 
be  derived  from  the  Hamiltonian  as  [29] : 


d/7  _  dH  ,  dH  J 

dH  —  >  — dqj+ >  -=r — dpt  H — z-—dt 

i  d([i  •  opt  ot 


(24) 


After  substituting  Equation  23  into  Equation  24, 


dH  =  ^4idpi  ~  Y, 


(25) 


and  from  Equation  22  it  follows  as, 


dL 


dH  =  Yd  idpi~YPidcli  ~  ~^dt 


(26) 


After  equating  each  term  in  Equation  26  to  each  term  in  Equation  24, 


dH 


dH 


di  —  3  5  Pi  —  3 

opi  dcfa 


(27) 


Equations  27  are  the  canonical  equations  of  Hamilton,  which  are  2n  number  of  first-order 
differential  equations. 


2.7  Canonical  Transformations 

Hamiltonian  formalism  doesn’t  introduce  an  efficient  calculation  tool  by  itself.  The 
freedom  in  choosing  the  physical  quantities  as  coordinates  and  momenta  is  beneficial  for 
solving  difficult  dynamical  problems.  Changing  coordinates  helps  to  reveal  much  about  any 
dynamical  system.  In  Hamiltonian  mechanics,  the  Hamiltonian  is  constant  if  it  doesn’t  ex- 
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plicitly  depend  on  time,  and  the  momenta  are  constant  if  the  associated  coordinates,  which 
are  called  ignorable  or  cyclic  coordinates,  are  missing  in  the  Hamiltonian.  However,  it  is 
necessary  to  start  from  the  beginning  when  the  change  of  coordinates  is  required.  More¬ 
over,  for  Hamiltonian  formalism,  the  coordinates  can  be  changed  to  anything  as  long  as 
they  describe  the  system,  but  the  momenta  are  dictated  by  the  Lagrangian,  see  Equation  22. 
One  method  for  changing  coordinates  without  starting  from  the  beginning  by  conserving 
the  structure  of  Hamiltonian  dynamics  is  canonical  transformation  [64]. 

Assume  Equation  28  represents  new  canonical  coordinates,  and  Equation  29  represents 
new  momenta: 

Qi  =  Qi{qi,Pi,t )  (28) 


Pi  =  Pi(qi,Pi,t) 


(29) 


Also  assume  a  new  Hamiltonian  K  —  K(Qi1Pi1t),  and  Hamilton’s  equations  as: 


Q i 


dK 


(30) 


The  old  variables  <r/,  and  /?/,  and  the  new  variables  Qj  and  /J,  are  supposed  to  describe  the 
same  dynamical  system.  The  Hamilton’s  principle  can  be  used  to  show  that  the  variation 
of  both  integrals  equal  to  zero  as: 

r  n 

8  /  i^Pi4i~H)dt  =  0,  (31) 

■'  i=i 


r  N 

8  C£PiQi~H)dt  =  0,  (32) 

J  7=1 

where  L  can  be  deduced  from  Legendre  transform  as  L  =  (Li=i  Pi4i  ~  H)-  Although  the 
variations  of  both  integrals  equal  to  zero,  the  integrands  aren’t  equal.  The  total  time  deriva- 
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tive  of  an  arbitrary  function  F  separates  two  integrands  from  each  other  as: 


s  I E Pi4i  ~  H ( Pi , q»t)  -  Y,piQi 

dF 

+K(Qi,Pi,t )  -  -7~]dt  =  0 


(33) 


The  variation  vanishes  at  end  times  as: 


5  f  2  ^ dt  =  S(F(t2)-F(t 0)  =  0,  (34) 

Jt\  dt 

where  F  is  the  generating  function,  which  is  a  function  of  independent  2n  variables.  Ta¬ 
ble  6  shows  the  four  possible  forms  of  the  generating  function  F.  The  generating  function 
F  needs  to  be  known  in  order  to  apply  the  transformation.  F  can  be  obtained  from  the 
relationship  between  old  and  new  coordinates,  and  Table  6,  see  Wiesel  [64]. 

Table  6.  Canonical  Transform  Relations 


Fi(q,Q,t) 

p2  (q,P,t) 

Fs{p,Q,t) 

Fa  (p,P,t) 

d  — 

Pl  ~  r)q, 

n  ■  —  M 

Pl  ~  dqi 

dF3 

qi~~w, 

_  _dF4 
h  ~  dpi 

p.  _  dF] 
p  ~  dQi 

n  —  <LFi 

-  dPi 

p  —  dFj, 

P  ~  dQj 

o  —  — 

y 1  _  dPj 

The  new  Hamiltonian  K  and  the  old  Hamiltonian  H  are  related  by  Equation  35  as: 

dF\ 

K(Q,P,t)=H(q,p,t)  +  -^  (35) 

2.8  Hamilton-Jacobi  Theory 

The  easier  solution  to  a  dynamical  system  can  be  accomplished  when  the  more  coor¬ 
dinates  are  missing  in  the  Hamiltonian.  Moreover,  any  momenta  missing  from  the  new 
Hamiltonian  dictates  that  its  conjugate  coordinate  is  constant.  If  the  new  momenta,  K, 
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equals  to  zero,  all  the  new  coordinates  and  momenta  are  constant  [64].  Equation  36  shows 
that  the  new  Hamiltonian  equals  to  zero,  which  is  called  the  Hamilton-Jacobi  equation,  and 
Equation  37  shows  that  the  new  coordinates  and  momenta  are  constant  [29]. 


H(qi,Pi,t)  + 


dF 

~dt 


0 


(36) 


Qi 


Pi 


dK 

dPi 

dK 


0 

0 


(37) 


It  is  favorable  to  use  a  generating  function  which  is  a  function  of  the  old  coordinates, 
qi,  and  the  new  momenta,  P,.  The  generating  function,  F2,  from  Table  6  is  such  a  function. 
Therefore,  Equation  36  becomes  [29]: 


H 


dF2  \  dF2  _ 
Qi:  — T  I  +  —  0 


dqi 


dt 


(38) 


Equation  38  is  known  to  be  Hamilton’s  principal  function.  The  solution  to  this  equation  is 
denoted  by  S  as  [29] : 

S  =  S(q\  •  •  ■  qn,  CC\  •  •  •  (39) 

where  a,  are  n  number  of  constants  of  the  integration  because  S  is  the  solution  to  the  first- 
order  differential  equation.  Therefore,  one  of  the  variables  of  the  solution  has  to  be  an 
arbitrary  constant  added  to  the  solution  S.  Using  the  analogy  of  the  physical  description 
of  the  generating  function,  the  constants  of  the  integration,  a,,  can  be  chosen  as  the  new 
momenta,  Pj  [29] : 

Pi  =  0Ci  (40) 
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After  applying  transformation  for  the  fo  generating  function  from  Table  6  [29] : 


dS(qi,(Xi,t) 

dqi 


Continuing  with  the  second  half  of  the  transformation  [29]: 


Qi  =  A 


dS(qj,  (Xj,t) 

c)al 


(41) 


(42) 


Both  a,  and  /3,  can  be  obtained  by  evaluating  partial  derivatives  with  the  known  initial 
conditions.  The  Hamilton-Jacobi  equation  can  be  solved  as  [29]: 


q  =  q((Xi,Pi,t) 


(43) 


The  Hamilton’s  principal  function  helps  to  transform  the  old  variables  to  new  constant 
coordinates  and  momenta.  A  solution  to  the  dynamical  system  is  obtained  when  solving  the 
Hamilton-Jacobi  equation.  The  relationship  between  the  Hamilton’s  principal  function  and 
the  Hamilton’s  principle  can  be  shown  by  examining  total  time  derivative  of  the  Hamilton’s 
principle  function,  S  [29] : 


dS 

dt 


v.  dS 
i  dqt 


ds 

~di 


(44) 


After  inserting  Equation  36  and  Equation  41,  Equation  44  becomes: 


dS 

dt 


Y^Pidi-H  =  L 


(45) 


Therefore,  a  constant  differs  the  Hamiltonian  principle  function  from  Hamilton’s  principle 
as: 

S  =  I  Ldt  +  constant  (46) 

For  Hamilton’s  principle,  the  solution  to  a  dynamical  system  is  obtained  by  solving  the 
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definite  integral  of  Lagrangian,  L.  Equation  46  represents  the  same  integral  in  indefinite 
form  which  can  be  used  to  solve  a  dynamical  system. 

Hamilton’s  principle  function  can  be  separated  in  two  parts  when  the  Hamiltonian 
doesn’t  explicitly  depend  on  time  as  [29] 


S(qt,  0Ci,t )  =  W (<?/, «/)  -  oqf,  (47) 

where  W(qi,  a,-)  is  the  first  part  of  the  solution,  which  depends  on  the  old  coordinates 
and  the  new  momenta,  and  appears  only  when  the  Hamiltonian  is  constant.  VE(g,-,cq)  is 
called  Hamilton’s  characteristic  function.  Equation  47  is  the  solution  to  the  Hamilton- 
Jacobi  equation  [29]: 

f+«te.f)  =  0  (48) 

After  substituting  Equation  47,  the  Equation  48  becomes  [29]: 

dW 

H{qi,  -5—)  =  Of!,  (49) 

where  one  of  the  constants  of  integration,  oq,  equals  to  the  constant  value  of  the  Hamilto¬ 
nian,  H.  The  new  Hamiltonian,  K,  equals  to  oq  because  H  doesn’t  depend  on  time  [29]: 

K=a  1  (50) 


After  calculating  Hamilton’s  equations  for  the  new  Hamiltonian,  K  [29]: 


Pi 


Qi 


dK 

dQi 

dH 

dai 


0,  Pj  =  at 

1  i  =  l 


=  0  1 


(51) 
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After  denoting  the  generating  function  by  W iq,.  Pt).  the  solution  becomes  [29]: 


Q  t 


Qi 


t  +  P  i  = 


dw 

da\ 


i  ^  1 


(52) 


where  Q\  is  not  a  constant  of  the  motion,  and  time  ,  t,  is  a  coordinate  and  the  Hamiltonian 
is  its  conjugate  momenta  [29]. 


2.9  Action-Angle  Variables 


The  periodic  motion  is  of  great  importance  in  physics.  The  frequencies  of  the  motion  is 
mostly  more  desirable  than  the  other  properties  of  the  orbit.  A  variation  of  the  Hamilton- 
Jacobi  process  is  used  to  obtain  the  frequency  of  the  periodic  motion.  The  action  variables, 
J,  are  chosen  instead  of  the  new  momenta  of  the  Hamilton- Jacobi  equations.  /,  are  inde¬ 
pendent  functions  of  the  a,  of  the  Section  2.8.  Figure  7  shows  an  example  of  the  two  types 
of  periodic  motion.  For  oscillation,  the  system  repeats  its  track  for  every  point  because  q 
and  p  return  to  its  original  values  after  one  period.  For  rotation,  the  position  coordinate  q 
is  an  unbounded  angle  of  rotation,  which  grows  indefinitely  with  a  period  of  qa,  whereas  q 
is  bounded  for  oscillation  [29] . 

Ji  can  be  defined  as  [29]: 

Ji  =  j>  pidqi  (53) 

After  applying  the  canonical  transformation,  Equation  53  becomes  [29]: 


Ji 


dW  (qj,  CCj) 
dqi 


dqi, 


(54) 


where  the  solution  to  the  Equation  54  is  Jj  as  independent  functions  of  a,.  Substituting  /, 
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Oscillation  Rotation 


Figure  7.  The  Orbit  of  Phase  Space  for  Periodic  Motion  [29] 

for  a,  yields  the  characteristic  function  W  as  [29]: 


W  =  W(qi"-qn;Ji"-Jn) 


Also,  the  angle  variables,  C0(,  are  the  conjugate  coordinates  to  the  J,  [29]: 


Pi  = 


co,  = 


dW 

dqt 

dW 

~dJi 


The  Hamiltonian  is  a  function  of  only  J,  because  w,-  are  cyclic  [29]: 


H  =  H(Jl--Jn)  =  al 


The  Hamilton  equations  of  motions  for  the  new  variables  become  [29] : 


jt  = 

CO,  = 


dH 

dojj 

dH 


(55) 


(56) 


(57) 


(58) 
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where  v,-  are  constant  functions  of  J,.  The  solution  becomes  [29]: 


Ji  =  constant 


(Oi  =  Vit  +  pi 


(59) 


where  the  solution  doesn’t  introduce  any  advantage  over  the  a,  coordinates  from  the  Sec¬ 
tion  2.8.  The  advantage  of  this  process  is  that  the  Action-Angle  variables  transformation 
yields  the  frequencies,  V/,  of  the  periodic  motion  without  finding  a  complete  solution.  If  a 
system  is  periodic,  the  energy  of  the  system  can  be  defined  in  terms  of  /,  by  Equation  53. 
Next,  the  frequencies,  v,,  of  the  system  can  be  obtained  by  Equation  58.  The  ft),  are  conve¬ 
niently  called  angle  variables  due  to  Equation  59  [29]. 

The  fact  that  the  frequency  v,-  are  related  to  </,  can  be  proved  by  investigating  the  change 
in  one  of  the  angle  variables,  ft),,  when  one  of  the  coordinates,  qt,  completes  its  one  cycle 
of  oscillation  or  rotation.  The  Equation  60  shows  that  the  change  in  the  angle  variable  is 
caused  by  an  incremental  increase  in  one  of  the  coordinates,  qj  [29]: 


A  ft),  = 


(60) 


where  5  ft),  is  the  infinitesimal  change  due  to  the  incremental  increase  in  one  of  the  coordi¬ 
nates,  qj  [29]: 

8(Oi—^r—^-dqj  (61) 

dqj 

Also,  by  Equation  56  and  Equation  61,  Aft),  can  be  written  as  [29]: 


A  ft), 


d(0 j 

3  dq  j  — 

dqj 


d2W  ^ 

d^JdJ,dqj 


(62) 


The  Jj  in  the  denominator  can  be  taken  out  of  the  integral  because  J,  are  independent  func- 
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tions  of  only  a,  [29]: 


Aw, 


d 

dJi 


By  Equation  53,  A  ft),  becomes  [29]: 


Aft), 


dJi 

dJi 


1 


i  =  j 


=  Sij  i  i 


If  t i  is  the  period  of  the  motion  related  to  qi.  Equation  65  becomes  [29]: 


(63) 


(64) 


A  ft)/  =  1  =  V/T/ 

,  (65) 

V/=  - 
T/ 

2.10  General  Canonical  Transformations 

It  is  not  an  easy  task  to  create  a  generating  function  for  every  canonical  transformations. 
A  different  approach  to  canonical  transformations  can  be  initiated  as  [64] : 


xT  =  {qu  Pi} 


(66) 


Next,  the  Hamilton’s  equations  of  motion  becomes  [64]: 


4i  = 


Pi  = 


c)H 

dpi 

dH 

dqi 


(67) 
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Equation  67  can  be  defined  in  terms  of  phase  space  state  vector  by  Equation  66  [64]: 


x  —  Z 


Z  = 


dH 

dx 

0  / 
-1  0 


(68) 


where  Z  satisfies  the  properties  [64]: 


-z  =  zl  =z 


(69) 


Nothing  new  has  been  introduced  except  representing  the  Hamilton’s  equations  in  terms 
of  the  phase  space  state  vector.  The  transformation  of  the  phase  space  from  the  old  vari¬ 
ables,  x,  to  the  new  variables,  y,  can  be  accomplished  as  [64]: 

X  =  f(y)  (70) 

Then,  the  new  Hamiltonian,  K ,  becomes  [64]: 


m = mf  w), 


(71) 


where  the  old  Hamiltonian,  H,  doesn’t  explicitly  depend  on  time.  The  transformation,  /,  is 
valid  as  long  as  it  is  a  canonical  transformation.  Therefore,  the  conditions  which  make  /  a 
canonical  transformation  should  be  examined.  The  new  coordinates  and  momenta  should 
comply  with  the  Hamilton’s  equations  as  [64]: 


y  =  Z 


dK 

dy 


(72) 


Converting  Equation  72  to  Equation  68  may  reveal  these  conditions.  Equation  73  shows 
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the  time  derivative  of  the  transformation,  and  Equation  74  represents  the  gradient  of  the  new 


Hamiltonian  [64]: 


(73) 


c)K  d  ..  fdfVdH 
lTy=TyH(m  )={i)  to 


(74) 


Combining  these  equations  with  Equation  72  yields  [64]: 


.  ,3/  l.  df  rdH 

y={Ty)  X  =  Z(Ty]  to 


(75) 


Then,  Equation  75  yields  i  as  [64]: 


•  tdf>7tSf>TSH 

X=(Ty)Z(Ty]  to 


(76) 


Comparing  Equation  68  to  Equation  76  immediately  yields  [64]: 


[dy)Z[dy)  Z’ 


(77) 


where  the  first-order  partial  derivative  is  called  the  Jacobian  matrix  7.  The  7  matrix  is 
symplectic  when  it  satisfies  [64]: 

7Z7r  =  Z  (78) 

Therefore,  the  transformation,  /,  is  a  canonical  transformation  if  its  Jacobian  matrix  is 
symplectic. 


2.11  Orbit  Perturbations  for  Near-Earth  Satellites 

Orbit  perturbations  are  the  deviations  from  the  two-body  orbit  motion.  Two-body  mo¬ 
tion  assumes  that  secondary  body  orbits  around  a  point  of  mass  or  spherically  symmetric 
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sphere  [15].  However,  Earth  is  neither  a  point  of  mass  nor  spherically  symmetric.  The 
equatorial  bulge,  continental  blocks,  ocean  basis,  and  mountain  ranges  result  in  gravita¬ 
tional  field  deviations  from  the  two-body  point  of  mass.  Another  important  perturbation 
for  near-Earth  satellites  is  air  drag.  Air  drag  effect  is  bigger  when  a  satellite  has  smaller 
mass  with  bigger  cross  section.  Air  drag  acceleration  is  a  non-conservative  force,  which 
determines  the  orbit  lifetime  of  a  satellite  in  Low  Earth  Orbit  (LEO)  [64].  Table  7  repre¬ 
sents  three  categories  for  Earth  orbits.  LEO  satellites  are  primarily  affected  by  air  drag  and 
non- sphericity  of  Earth,  and  both  MEO  and  GEO  satellites  are  mainly  perturbed  by  solar- 
radiation  pressure  and  third-body  effects  [60].  This  work  includes  air  drag,  and  sectoral  and 
tesseral  harmonics  perturbation  accelerations.  The  effects  of  zonal  harmonics  perturbation 
are  included  in  the  periodic  orbit,  see  Section  3.2. 

Table  7.  The  Three  Categories  of  Earth  Orbits  [60] 


Orbit 

Altitude  (km) 

Low  Earth  Orbit  (LEO) 

h  <  800 

Medium  Earth  Orbit  (MEO) 

800  <  h  <  30,000 

Geostationary  Orbit  (GEO) 

h  =  35,780 

The  gravitational  field  equation,  which  is  also  called  Poisson’s  equation,  is  the  funda¬ 
mental  partial  differential  equations  for  gravitational  fields  [64]: 


V2V  =AnGp,  (79) 

where  G  is  the  gravitational  constant,  p  is  mass  distribution,  V2  is  Laplacian  operator,  and 
V  is  gravitational  potential  function.  Equation  79  defines  the  gravitational  potential  of  a 
sphere.  The  gravitational  potential  can  be  calculated  by  Equation  79  when  mass  distri¬ 
bution,  p,  is  provided.  However,  the  gravitational  potential  outside  the  sphere  is  needed. 
Equation  80  represents  the  infinite  series  of  the  gravitational  potential  function  outside  a 
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Zonal  Harmonics 


2,0  3,0  4,0  5,0  n,m 


Sectorial  Harmonics 


2,2  3,3  5,5  n'm 


Tesseral  Harmonics 


3,1  3,2  4,1  4,2  n,m 


Figure  8.  The  Visual  Depictions  of  Spherical  Harmonics  [60] 


sphere,  which  has  a  radius  of  P©  [64]. 

vm,<m  =  - 7 £  £  (^r“C(«*e) 

n=0m=0  ®  (80) 

x  ( Cnmcosmty  +  Snmsinmty) 

where  /l  is  the  gravitational  parameter,  n  and  m  are  the  order  and  degree  of  the  expansion, 
respectively,  7?©  is  the  radius  of  Earth,  P™  are  the  associated  Legendre  functions,  the  0 
is  the  geocentric  latitude,  0  is  the  longitude,  and  Cnm  and  Snm  are  gravitational  constants, 
which  can  be  obtained  from  the  gravitation  models.  Both  Cnm  and  Snm  define  the  shape  of 
the  gravitational  field  [64].  All  spherical  harmonics  in  the  gravitational  potential  field  of 
Earth,  which  are  zonal,  tesseral,  and  sectoral  harmonics,  can  be  represented  by  Equation  80 
[60].  Figure  8  shows  the  visual  depictions  of  three  types  of  spherical  harmonics. 

Newtonian  point  of  mass  potential  can  be  obtained  from  Equation  80  when  m,n  =  0. 


39 


The  first  term  is  free  of  longitude  and  latitude  terms  because  cos  m(j>  =  0,  sin  nW  =  0, 
Coo  =  1,  and  the  associate  Legendre  function  lJ[]  =  1 . 

V  =  --,  (81) 

r 

where  V  is  the  Newtonian  point  of  mass  potential. 

The  geopotential  expansion  yields  zonal  harmonics  when  n  is  the  order  and  m  —  0 
is  the  degree.  The  smallest  zonal  harmonic,  which  has  n  =  1  order  and  m  —  0  degree, 
is  considered  to  be  zero  because  it  shifts  the  center  of  mass  of  Earth  parallel  to  North- 
South  axis.  Any  frame  which  doesn’t  have  its  origin  at  the  center  of  mass  is  not  practical. 
Therefore,  the  first  zonal  harmonic  equals  to  zero,  J\  =  0,  when  the  origin  is  defined  at  the 
center  of  Earth.  The  second  zonal  harmonic  specifies  the  oblateness  of  Earth  at  the  equator. 
For  the  Earth’s  potential  expansion,  the  J2  is  second  to  the  Newtonian  point  mass  potential 
in  strength.  The  Newtonian  potential  is  103  times  bigger  than  the  the  72.  The  contribution 
of  J2  term  to  the  gravitational  potential  can  be  represented  as  [64] : 

V20  =  ^S^(3cos2d-l),  (82) 

where  J2  —  —C2o  —  0.001082.  For  the  second  zonal  harmonic,  J2,  mass  concentration  is 
bigger  at  the  equator  than  the  poles.  The  higher  order  zonal  harmonics  add  more  terms  to 
the  potential,  which  model  the  irregularities  of  the  mass  distribution  along  the  longitudes. 
See  Figure  8  for  the  visual  depiction  of  zonal  harmonics  [64]. 

The  sectoral  harmonics  are  the  terms  with  equal  orders  and  degrees  in  the  potential 
expansion,  n  =  m.  They  are  functions  of  only  longitude,  and  they  change  signs  across 
the  longitudes.  The  first  sectoral  harmonic,  n  =  m  —  1,  changes  the  location  of  center  of 
mass  away  from  the  origin  of  a  reference  frame  to  another  one  in  the  equatorial  plane. 
Therefore,  Sn  =  Cn  =0  because  it  is  always  desirable  to  have  center  of  mass  at  the  origin 
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of  a  reference  frame.  The  second  sectoral  harmonic,  n  =  m  =  2,  contributes  to  the  potential 
as  [64]: 

V22  «  C22COS  2 (j)  +  S22sin  2 (j) ,  (83) 

where  (j)  is  the  longitude.  The  second  sectoral  harmonics,  C22  and  S22,  are  on  the  order  of 
10  6 .  J2  is  103  times  bigger  than  C22  and  .So 2-  The  mass  concentrations  of  Earth  at  Eurasia 
and  the  Americas  are  bigger  than  Pacific  and  Atlantic  ocean  basis  due  to  the  second  sectoral 
harmonics.  See  Figure  8  for  the  visual  depiction  of  sectoral  harmonics  [64]. 

The  tesseral  harmonics  are  the  case  when  n^m^  0  in  the  geopotential  expansion.  They 
create  a  sphere  with  rectangular  square-tiled  boards  on  it.  The  number  of  bands  in  latitude 
is  n  —  m+  1  because  of  P"'(cos  Q)  dependence.  The  terms  ,  Cnmcos  m(j)  and  Snmsin  m(j), 
disappear  for  2m  meridians.  Therefore,  there  are  2m  bands  in  longitude.  See  Figure  8  for 
the  visual  depiction  of  tesseral  harmonics  [64,  60]. 

The  gravity  models,  such  as  EGM96,  include  the  coefficients,  Cnm  and  Snm.  The  gravity 
models  are  built  from  terrestrial  or  satellite  measurements.  The  satellite  measurements  are 
orbit  dependent.  Therefore,  the  accuracy  of  a  gravity  model  depends  on  the  number  of 
satellites  with  different  orbits.  This  effort  uses  the  EGM96  gravity  model.  The  EGM96  is 
developed  by  the  University  of  Texas  at  Austin,  the  Defence  Mapping  Agency,  Ohio  State 
University,  and  the  Goddard  Spaceflight  Center.  The  EGM96  includes  30  different  satellite 
measurements  and  terrestrial  gravity  measurements.  It  is  a  complete  model  with  degree  and 
order  360  [60].  However,  this  work  uses  only  order  and  degree  20  because  it  provides  high 
accuracy  with  low  computational  power  requirement.  Figure  9  represents  the  variations  in 
the  geopotential. 

The  air  drag  is  the  other  perturbation  included  by  this  work.  The  air  molecules  behave 
individually,  which  is  called  the  free  molecular  flow  regime,  at  the  altitude  of  satellites. 
The  air  molecules  impact  the  motion  of  the  satellites.  Equation  84  represents  the  drag  law 
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Europe 
&  Africa 


Figure  9.  GRACE  Gravity  Model  [53] 


Temperature  Variations  in  the  Atmosphere 


Figure  10.  Temperature  Variations  in  the  Atmosphere 


[64]: 

ad  =  -\— pVV,  (84) 

2  m 

where  Cd  is  the  drag  coefficient,  which  is  2.2  for  flat  plane  models,  and  between  2.0  and 
2. 1  for  spheres  for  the  most  part,  p  is  the  atmospheric  density,  A  is  the  cross-sectional  area 
of  the  satellite,  V  is  the  velocity  of  the  satellite  relative  to  the  air  molecules,  and  m  is  the 
mass  of  the  satellite.  The  drag  coefficient,  Cd,  specifies  the  vulnerability  of  a  satellite  to 
the  air  drag,  and  it  is  dimensionless.  The  atmospheric  density  at  the  altitude  of  a  satellite 
is  represented  by  p,  which  is  difficult  to  calculate  due  to  the  interaction  between  Earth’s 
upper  atmosphere  and  Sun.  The  cross-sectional  area,  A,  is  also  difficult  to  calculate  because 
the  attitude  of  the  satellite  must  be  known.  Moreover,  it  is  almost  impossible  to  determine 
the  attitude  of  a  tumbling  satellite  [60] . 

Solar  flux  and  geomagnetic  storms  heat  up  the  thermosphere.  Therefore,  the  upper  at- 
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mosphere  expands,  which  increases  the  number  of  air  molecules.  During  Solar  maximum 
the  emission  of  Extreme  Ultraviolet  radiation  (EUV)  increases  because  Solar  flare  and 
Coronal  Mass  Ejection  (CME)  events  increase  during  Solar  maximum.  Charged  particles 
coming  from  the  night  side  of  the  Earth’s  atmosphere  during  geomagnetic  storm  interact 
with  air  molecules,  and  increase  their  energy,  which  expands  the  upper  atmosphere.  Ge¬ 
omagnetic  storms  cause  delayed  expansion  in  the  atmosphere  because  storms  first  hit  the 
day  side  of  the  Earth’s  atmosphere.  Both  Solar  flux  and  geomagnetic  storms  increase  the 
atmospheric  density  in  the  upper  atmosphere.  Figure  10,  which  is  plotted  using  Mass  Spec¬ 
trometer  -  Incoherent  Scatter  (MSIS-E-90)  model  data,  shows  the  temperature  change  in 
the  atmosphere  due  to  the  diurnal  variations  and  the  Solar  cycle,  and  Figure  11,  which  is 
plotted  using  SGP4  and  TLE  prediction  data,  represents  the  impact  of  air  drag  on  the  Delta 
1  rocket  body  over  138  days,  which  has  an  altitude  of  482  km. 
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Figure  11.  Air  Drag  Effect  on  the  Delta  1  Rocket  Body 

Although  the  atmospheric  density  is  unpredictable,  it  is  very  important  for  the  accuracy 
of  the  orbit  predictions.  Many  atmospheric  models,  which  are  either  static  or  time-varying, 
has  been  developed.  Every  model  is  either  based  on  physical  models,  which  are  built 
by  combining  conservation  laws  and  atmospheric-constituent  models  or  developed  from 
in-situ  measurements  and  satellite-tracking  data.  Every  atmospheric  model  is  different  in 
terms  of  speed,  accuracy,  and  applicability.  Therefore,  there  is  no  model  which  provides 
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best  results  for  all  applications.  This  effort  uses  the  1962  Standard  Atmosphere.  The 
routine  for  the  atmosphere  model  requires  the  altitude  of  each  layer,  molecular  weight,  and 
molecular  temperature.  Moreover,  the  temperature  within  each  layer  changes  linearly,  see 
Regan  and  Anandarskarian  [27].  It  is  an  ideal,  static  atmospheric  model  at  a  latitude  of  45° 
during  moderate  Solar  activity  [60]. 


2.12  Numerical  Integration  Methods 


There  are  many  numerical  integration  methods  which  have  been  developed  to  solve 
ordinary  differential  equations.  Most  of  these  methods  have  been  successfully  applied  in 
celestial  mechanics.  However,  there  is  no  numerical  integration  method  which  yields  the 
best  solution  to  every  problem  pertaining  to  the  motion  of  satellites.  The  most  important 
numerical  integration  methods  for  orbit  computations  are  Runge-Kutta  methods,  multistep 
methods,  and  extrapolation  methods.  Runge-Kutta  methods  are  easy  and  applicable  to  wide 
range  of  different  problems.  Multistep  methods  are  very  efficient,  but  they  need  to  store 
previous  data  points.  Extrapolation  methods  are  highly  accurate  [56]. 

This  effort  uses  4?/'-order  Hamming  predictor-corrector  method,  which  is  basically  a 
multistep  numerical  integrator.  It  requires  four  values  of  the  state  vector.  However,  only 
the  initial  conditions  are  know  at  the  beginning.  The  other  three  values  can  be  obtained  by 
a  process  called  Picard  iteration.  Equation  85  shows  the  Picard  iterative  process: 


yn+i{x)=yn+  f(t,yn(t))dt 

Jx°  (85) 

yn  0)  =yn 


If  the  Picard  iteration  converges,  the  four  values  of  the  state  vector  are  extrapolated  to 
give  the  next  value  for  the  state  vector,  which  is  the  predictor  part.  Next,  the  extrapolated 
state  is  corrected  using  better  values  of  the  equations  of  motion,  which  is  the  corrector 
part.  Although  predictor-corrector  methods  require  complex  algorithms,  they  are  favorable 
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because  of  their  stability. 


2.13  KAM  theory 

The  French  mathematician  and  physicist  Henri  Poincare  found  out  that  3BP  is  unsolv- 
able  because  of  small  divisors  problem.  He  suggested  that  N-Body  Problem  is  mathemati¬ 
cally  unsolvable.  Birkhoff  and  Smale  proved  Poincare’s  earlier  findings,  and  they  suggested 
that  some  nonlinear  systems  aren’t  solvable.  The  primary  question  was:  When  does  a  sys¬ 
tem  behave  chaotically?  The  KAM  theory,  which  was  announced  by  A.N.  Kolmogorov  in 
1954,  and  proved  by  J.  Moser  and  V.I.  Amol’d  in  the  early  1960’s,  is  the  solution  to  this 
question  from  the  Hamiltonian  aspect.  The  KAM  theory  explains  what  happens  when  an 
integrable  Hamiltonian  is  perturbed.  It  is  developed  to  overcome  the  difficulties  in  pertur¬ 
bation  theory  about  small  divisors  [50].  The  KAM  theory  states  that  an  integrable,  which 
means  n  constants  of  the  motion  are  known,  Hamiltonian  has  a  phase  space  motion  which 
lies  on  a  ^-dimensional  torus  in  2/7-dimensional  phase  space,  where  n  is  the  number  of  in¬ 
dependent  coordinates  [47,  6].  The  qua  si-periodic  motion  can  be  represented  by  action  and 
angle  variables,  see  Section  2.9.  The  quasi-periodic  orbits  describe  integrable  motion  on 
the  invariant  torus.  If  any  trajectory  is  on  an  invariant  torus  in  phase  space,  it  will  stay  on 
the  torus  [5].  The  quasi-periodic  motions  have  n  number  of  frequencies,  which  is  dictated 
by  the  Hamilton-Jacobi  theory  [65]. 

The  fundamental  equation  for  the  KAM  theory  can  be  shown  as  : 

H(l:G)=H0(l)  +  eHl(l:d),  (86) 

where  H  is  the  perturbed  Hamiltonian,  H\  is  the  perturbing  Hamiltonian,  H0  is  the  in¬ 
tegrable  Hamiltonian,  e  is  a  small  real  number,  e  <<  1,  and  I  and  0  are  action-angle 
variables.  H0  and  H\  must  be  real  analytic  functions,  which  are  infinitely  differentiable 
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Figure  12.  A  Visual  Depiction  of  1-Torus  [50] 


functions.  The  action  variable  1  is  constant  and  0  is  a  linear  function  of  time  because  Ha 
is  only  a  function  of  I,  see  Equation  27.  Figure  12  represents  1-torus,  and  Figure  13  shows 
2-torus  where  the  action  variables,  I,  are  constant  and  the  angle  variables,  0,  are  functions 
of  time. 

If  the  solutions  for  the  perturbed  Hamiltonian  Equation  86  lie  on  ^-dimensional  tori, 
there  are  new  action-angle  variables  as  [57]: 


(87) 


Equation  87  represents  the  generating  function,  S,  for  the  Hamilton- Jacobi  transformation, 
see  Section  2.8  [57]: 


dS(l') 


<90 

0W) 


(88) 


New  Hamiltonian  is  a  function  of  only  I’  because  the  0  variables  are  cyclic,  and  the 
Hamilton- Jacobi  equation  becomes  [57]: 


(89) 
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Figure  13.  A  Visual  Depiction  of  2-Torus 

Power  series  expansion  in  £  can  yield  a  solution  to  the  generating  function,  S,  as  [57]: 


S  —  Sa  +  eS\  +  £2S2  H -  (90) 


Substituting  Equation  89  in  Equation  88  assuming  S0  —  I'-O  because  Equation  87  yields 
1  =  1'  and  0  =  0'  after  the  assumption  [57]: 


+  (I  +  £ 


<9*52 

do 

dS  i 

~d0 


+  •••) 

+  ...,0)=H\ V) 


(91) 


The  expansion  of  Equation  91  for  small  £  keeping  only  the  first-order  terms  yields  [57]: 


H0(  l')+£ 


dHc  dS i 
~W~dG 


+  £Hi(t  ,0)  =H'(t) 


(92) 


Then,  H i  (I',  0)  and  Si  (I;,  0)  can  be  represented  as  Fourier  series  [57]: 


Hi  =  ■  0) 

m  (93) 

Si  =  ^Si,m(I/)erp(/m  ■  0) 

m 

where  m  is  a  vector  with  n  rows  consisting  of  integers.  Substituting  Equation  93  in  Equa- 
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Figure  14.  Elliptic  Islands  and  Hyperbolic  Fixed  Points 


tion  92  yields  [57]: 


(94) 


where  va  is  the  unperturbed  ^-dimensional  frequency  vector  by  Equation  59.  The  infinite 
sum  must  converge  to  provide  a  solution.  Moreover,  the  same  condition  must  be  satisfied 
for  S2,  S3,  •••  if  £  is  expanded  to  higher  orders.  Although  a  method  of  successive  ap¬ 
proximations  to  S  has  been  shown  here,  the  proof  of  KAM  theory  depends  on  complex 
successive  approximations  which  converge  much  faster. 

The  convergence  of  the  sum  depends  on  the  denominator,  m  ■  V0(I/).  This  problem  is 
called  the  small  divisors  problem.  The  action  variables,  I,  specify  the  resonant  tori,  which 
behave  chaotically  when  £  >  0,  for  the  unperturbed  system  when  m  ■  v0(I')  =  0.  On  the 
other  hand,  the  nonresonant  tori  must  satisfy  the  Equation  94  [57]: 


m  ■  v|  >  Ar(v)|m|  91+1) 


(95) 


where  m  are  vectors  consisting  of  integers,  the  absolute  value  of  m,  |m|,  is  the  sum  of  all  m 
vectors  up  to  number  n,  and  K(v)  >  0  is  an  arbitrary  number  independent  of  m.  Moreover, 
m  can’t  be  the  zero  vector  for  nonresonant  tori  [57].  Therefore,  it  is  highly  likely  that  a 
perturbed,  nearly  integrable,  and  periodic  system  is  described  by  an  invariant  tori  in  phase 
space  [1]. 

Another  way  to  describe  the  resonant  tori,  which  helps  to  visualize  it,  is  by  winding 
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number,  a.  Figure  13  can  help  to  visualize  of  the  ratio  of  two  frequencies,  Vi  and  V2. 
Assume  Vi  is  the  frequency  seen  from  top  view,  and  V2  is  the  frequency  seen  from  cross 
sectional  view.  If  the  winding  number,  a  =  V1/V2,  is  a  rational  number,  the  tori  associated 
with  this  winding  number  are  called  resonant  tori.  The  tori  with  irrational  winding  number, 
a ,  are  called  nonresonant  tori.  For  nonresonant  tori,  “almost  all”  orbits  are  preserved 
because  “the  KAM  theory  restore  a  measure  of  continuity  to  chaos”  [50].  Figure  14,  which 
consists  of  two  surface  of  section  plots,  represents  what  happens  if  a  resonant  torus  breaks 
up  when  a  perturbation  is  introduced  to  the  system.  According  to  the  Birkhoff’s  theorem, 
alternating  elliptic  and  hyperbolic  fixed  points  are  formed  when  the  resonant  torus  breaks 
down.  The  two  stable  elliptic  islands,  which  is  represented  by  O ,  and  the  two  unstable 
hyperbolic  points,  which  is  represented  by  X,  are  shown  in  Figure  14. 

2.14  KAM  Theory  Applications 

After  KAM  theory  was  published,  astronomers  applied  it  to  astronomical  models  be¬ 
cause  the  motions  in  the  solar  system  are  comparatively  bounded.  However,  application 
of  KAM  theory  to  astronomical  models  didn’t  yield  promising  results  because  of  the  lim¬ 
itation  on  the  size  of  the  perturbation  parameter,  e,  which  is  the  ratio  of  masses.  In  1963, 
Arnold  endeavored  to  establish  the  existence  of  KAM  tori  for  NBP,  which  he  was  partially 
successful.  The  primary  question  was  whether  there  were  the  initial  position  and  velocities 
of  the  bodies  that  keep  the  distance  of  the  bodies  from  each  other  bounded  for  all  the  time 
in  the  NBP  [2].  A  combination  of  KAM  theory  and  computer-assisted  techniques  named 
interval  arithmetic,  which  is  a  computational  technique  that  controls  the  errors  introduced 
by  numerical  computations  in  a  special  way,  was  successfully  used  to  prove  the  existence 
of  KAM  tori  for  the  Restricted,  Planar,  Circular  3-Body  Problem  (RPC3BP)  by  Celletti  and 
others.  In  1997,  Celletti  and  Chierchia  proved  the  existence  of  quasi-periodic  tori  with  a 
frequency  close  to  the  average  frequency  of  Ceres  for  the  Sun-Jupiter-Ceres  problem  [12]. 
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Locatelli  and  Giorgilli  proved  the  existence  of  KAM  tori  describing  the  secular  motions  of 
Jupiter  and  Saturn  for  the  obserbed  values  of  the  parameters  [48].  Celletti  and  Chierchia 
investigated  a  truncated  RPC3BP  model  for  Sun,  Jupiter,  and  Asteroid  12  Victoria.  The  in¬ 
variant  tori  bounding  the  motion  of  Victoria  was  built  successfully  for  the  astronomic  value 
of  the  Jupiter-Sun  mass  ratio  [13].  For  further  information,  see  a  brief  history  of  KAM  tori 
for  NBP  [14]. 

The  numerical  applications  of  KAM  theory  originated  from  Binney  and  Spergel’s  pa¬ 
per  which  spectra  obtained  from  Fourier  series  of  the  coordinates  of  numerically  integrated 
orbits  under  the  influence  of  galactic  potentials  can  be  expressed  as  sums  of  integer  mul¬ 
tiple  of  fundamental  frequencies.  Binney  and  others  realized  that  orbits  of  A-dimcnsional 
galaxies  defined  in  phase  space  can  be  represented  as  N- tori  [4,  6].  McGill  and  Binney  de¬ 
veloped  a  method  to  build  tori  for  Hamiltonians  including  general  gravitational  potentials. 
This  method  is  based  on  the  idea  of  distorting  the  analytic  tori  of  a  toy  Hamiltonian  into  the 
desired  tori  using  generating  functions  [18,  6].  Binney  and  Kumar  generalized  the  method 
for  obtaining  the  frequencies  and  angular  variables  related  to  the  tori  least-squares  fitted  to 
any  Hamiltonian  [3,  6].  Kaasalainen  and  Binney  further  refined  the  torus-fitting  process  of 
the  method  in  order  to  include  the  case  of  non-rotating  bar,  which  admits  two  major  orbit 
families  instead  of  one  [36,  6].  Kaasalainen  and  Binney  introduced  point  transformations 
into  the  method  in  order  to  solve  the  problem  of  a  toy  potential  being  too  different  from 
its  target  potential  [35,  6].  Then,  Kaasalainen  showed  that  the  method  can  be  applied  to 
globally  chaotic  regions  as  well  [34,  6] . 

The  only  work  related  to  this  effort  in  terms  of  the  application  of  KAM  theory  for  ob¬ 
jects  orbiting  Earth  has  been  done  by  Wiesel  and  his  masters  and  PhD  students.  Their 
efforts  are  represented  in  a  chronological  order.  Wiesel  showed  that  Earth  satellite  orbits 
under  the  full  geopotential  effect  are  likely  to  be  KAM  tori  in  a  reference  frame  rotating 
with  Earth.  First,  he  determined  the  frequencies  of  the  orbit  using  a  Laskar  frequency 
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algorithm.  Next,  he  obtained  Fourier  series  coefficients  from  numerically  integrated  tra¬ 
jectories.  Then,  He  fitted  numerically  integrated  orbit  to  the  multiple  Fourier  series,  which 
defines  the  torus,  using  least  squares  fitting  technique.  The  spectral  analysis  of  almost 
all  orbits,  excluding  chaotic  ones,  yielded  three  basis  frequencies.  Moreover,  Wiesel  con¬ 
cluded  that  other  dynamical  effects  can  be  added  as  perturbations  because  of  the  single 
point  construction  of  the  Hamilton-Jacobi  separation  of  variables  solution  [65].  His  next 
paper  related  to  the  application  of  KAM  torus  theory  to  Earth  satellites  is  the  one  that 
compares  two  different  KAM  torus  construction  algorithms,  which  are  the  least  squares 
fitting  of  a  KAM  torus  to  a  numerically  integrated  orbit  and  the  extraction  of  the  Fourier 
coefficients  of  the  KAM  torus  from  numerically  integrated  Fourier  transform.  An  efficient 
KAM  torus  construction  method,  which  yields  the  torus  that  passes  through  the  given  initial 
conditions,  is  important  for  applications,  such  as  formation  flight  of  satellites,  stationkeep¬ 
ing  satellite  constellations,  and  compressing  ephemeris  data  for  navigational  satellites  [66]. 
Derbis  tried  to  apply  the  KAM  theory  to  precise  satellite  data  from  the  GPS  satellites.  She 
applied  a  Fast  Fourier  Transform  (FFT)  to  obtain  and  identify  discrete  frequencies  from 
orbit  data.  However,  some  orbits  showed  inconsistencies  related  to  their  basis  frequencies. 
She  concluded  that  the  inconsistencies  in  the  third  frequency,  (03 ,  which  is  the  apsidal  re¬ 
gression  rate,  caused  by  the  resonance  in  orbits  of  the  GPS  satellites  [22].  Little  tried  to 
build  KAM  tori  from  orbit  data  of  the  Gravity  Recovery  Climate  Experiment  (GRACE) 
and  Jason- 1  satellites.  He  applied  a  modified  Laskar  Fourier  transform  algorithm  to  obtain 
the  basis  frequencies  of  the  orbits.  Although  the  residuals  of  the  KAM  torus  for  Jason- 1 
satellite  were  as  low  as  1  km  over  a  30  day  period,  the  KAM  torus  construction  process 
failed  for  GRACE  satellite.  Little  concluded  that  air  drag  and  errors  related  to  Fourier 
coefficients  resulted  in  the  poor  results  for  GRACE  satellite.  Craft  investigated  the  appli¬ 
cability  of  the  KAM  theory  for  the  satellite  formation  flight  in  the  full  geopotential  with  an 
order  and  degree  of  20.  The  KAM  torus  method  yielded  promising  results  for  the  satellite 


51 


formation  flight  applications.  Orbits  of  satellite  clusters  provided  with  drift  rates  in  the 
nanometer  to  micrometer  range  over  60  day  integration  interval  when  they  were  separated 
by  10  to  100  meters  [20].  Bordner  tried  to  apply  the  KAM  theory  to  the  GPS  satellite 
orbits.  His  research  questions  were  whether  the  KAM  theory  can  increase  the  accuracy 
of  GPS  satellites  and  whether  it  can  reduce  the  burden  on  GPS  operations.  However,  the 
methods  for  constructing  KAM  tori  failed  to  yield  desired  accuracy  for  operational  GPS  or¬ 
bits.  He  suggested  that  improved  methods  are  needed  to  deal  with  complexities  related  to 
the  maneuvering  operational  satellites  [6].  Wiesel’s  third  paper  related  to  the  KAM  theory 
compared  the  KAM  tori  built  from  2BP,  SGP4  model,  and  numerically  integrated  orbit  with 
a  degree  and  order  of  20.  He  showed  that  a  KAM  torus  is  similar  to  a  full  analytic  perturba¬ 
tion  theory  in  many  ways.  However,  for  a  KAM  torus,  frequencies  aren’t  approximations 
of  the  perturbation  series  expansions,  and  the  torus  is  built  numerically.  Moreover,  the 
comparison  of  numerical  integration  and  the  KAM  torus,  which  yielded  residuals  smaller 
than  4  m  over  10  years,  showed  that  trajectories  in  the  full  geopotential  are  KAM  tori  [68]. 
Yates  investigated  the  motion  near  a  reference  torus  in  order  to  compensate  for  errors  in  the 
reference  KAM  torus  due  to  dissipative  forces,  such  as  air  drag,  lunar/solar  mass  effects. 
He  suggested  that  routine  phase  angle  updates  and  stochastic  offsets  to  the  reference  torus 
can  improve  the  accuracy  of  the  KAM  torus  perturbed  by  dissipative  forces.  He  showed 
that  most  stochastic  effects  can  be  modeled  to  predict  the  motion  near  a  reference  torus. 
Moreover,  the  low  eccentricity  of  International  Space  Station  (ISS)  is  proved  to  be  a  big 
hurdle  in  estimating  the  torus  coordinate  offsets  [70].  Hagen  studied  the  effects  of  air  drag 
and  lunar  mass  perturbations  on  a  reference  KAM  torus.  He  showed  that  the  KAM  theory 
is  also  applicable  when  air  drag  and  lunar  mass  perturbations  are  included  in  the  system. 
However,  it  is  proved  that  an  accurate  Jacobian,  which  performs  the  linear  transformation, 
is  of  great  importance,  and  Jacobian  is  only  valid  when  the  reference  torus  coordinates 
are  tightly  aligned  with  the  perturbed  torus  coordinates  [30].  This  effort  includes  air  drag 
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as  perturbations  to  the  reference  KAM  torus.  Frey  studied  feasibility  of  constructing  the 
KAM  torus  using  TLE  and  SGP4.  Frey’s  work,  in  a  sense,  proved  that  the  KAM  torus 
can  be  constructed  from  observational  data.  Two  of  his  test  cases,  which  are  the  two  Delta 
rocket  bodies,  failed  due  to  air  drag  and  inaccuracies  of  the  TLE  data.  The  other  two  test 
cases,  which  are  the  Hubble  Space  Telescope  (HST)  and  the  Thor  rocket  body,  showed 
that  it  is  possible  to  extract  accurate  KAM  torus  basis  frequencies  although  they  suffered 
from  their  small  eccentricities.  Frey  concluded  that  small  eccentricities  can’t  be  modeled 
with  a  modified  Laskar  frequency  algorithm  developed  by  Wiesel  [28].  This  effort  uses  a 
different  method  which  builds  the  reference  KAM  torus  for  orbits  with  low  eccentricities 
using  periodic  orbits  and  Floquet  theory. 

This  effort  is  based  on  Wiesel’s  recent  paper,  “A  Theory  of  Low  Eccentricity  Earth 
Satellite  Motion”,  and  it  builds  on  the  results  of  previous  efforts  related  to  the  application 
of  KAM  theory  to  Earth  orbiting  satellites.  Wiesel  considered  the  Earth  satellite  motion 
in  terms  of  periodic  orbits  and  Floquet  theory.  Periodic  orbits  in  the  zonal  potential  are 
known  to  be  nearly  circular,  excluding  the  ones  with  critical  inclination  [58].  The  main 
advantage  of  applying  this  method  is  that  perturbations  are  in  the  order  of  10  5  instead  of 
10  3  because  periodic  orbits  includes  the  effects  of  the  Earth’s  oblateness.  The  perturba¬ 
tions  caused  by  sectoral  and  tesseral  potential  terms,  air  drag,  and  third  body  mass  effects 
can  be  added  to  the  periodic  orbit.  The  current  method  is  a  hybrid  of  numerical  and  an¬ 
alytic  methods.  It  has  numerical  sets  of  algorithms  that  may  match  numerical  integration 
in  their  accuracy.  Moreover,  Once  the  theory  files  are  constructed,  which  takes  approxi¬ 
mately  1  minute  per  orbit,  it  has  the  the  usual  advantages  of  general  perturbations  because 
it  provides  the  results  at  the  time  of  interest  without  having  to  perform  a  long  propagation 
[69], 
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2.15  Summary 


Space  surveillance  systems  are  under  great  stress  due  to  the  number  of  objects  orbiting 
the  Earth.  Current  launch  rates  and  operational  lifetimes  of  current  satellites  show  that 
space  surveillance  systems  need  to  cope  with  more  objects  soon.  The  increase  in  the  spatial 
density  of  debris  leads  to  an  increase  in  the  probability  of  collision.  Figure  15  represents 
the  growth  of  the  cataloged  satellite  population  during  the  past  15  years.  Therefore,  a  new 
method  which  is  as  accurate  as  numerical  methods  and  as  fast  as  analytical  methods  is 
needed.  The  KAM  theory  for  Earth  orbiting  objects  yield  promising  results  specifically 
for  space  debris  because  they  are  non-maneuverable,and  lightly  perturbed.  Once  a  KAM 
torus  is  built  for  a  non-maneuverable  object,  its  motion  will  be  bounded  by  the  torus  until 
a  dissipative  force  affects  its  motion.  Formation  flight  of  satellites  is  another  potential 
application  of  the  KAM  tori.  Craft  showed  that  tight  formations  yield  secular  drift  rates 
between  satellites  from  4  nanometers  to  1  micrometer  per  second  [20].  Moreover,  The 
KAM  tori  can  reduce  the  operational  burdens  of  GPS  because  the  ephemerides  provided  by 
the  orbital  tori  will  have  a  much  longer  useful  life  [6].  Indeed,  the  KAM  tori  has  numerous 
potential  applications  because  they  provide  accurate  predictions  with  long  lasting  validity 
for  the  time  of  interest  as  fast  as  an  analytic  method  can  provide.  However,  computationally 
expensive  long  integrations  for  building  a  torus,  and  perturbations  that  shift  a  reference 
torus  to  the  nearby  torus  have  been  the  main  problems  for  the  theory.  Hopefully,  this 
effort  mostly  overcomes  these  issues  and  the  KAM  torus  orbit  determination  will  serve  as 
a  prelude  to  converting  the  TLE  catalog. 
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Altitude  (km) 


Figure  15.  The  Growth  of  the  Cataloged  Satellite  Population  over  15  Years  [55] 
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III.  Methodology 


This  chapter  describes  the  data,  the  method,  and  the  analytic  process  used  to  perform 
this  study.  This  effort  is  intended  for  small  eccentric  orbits,  which  were  one  of  the  main 
issues  for  the  KAM  torus  construction  process  for  the  previous  efforts.  Although  this  effort 
has  been  built  on  the  previous  works,  it  introduces  a  completely  new  method  for  construct¬ 
ing  the  KAM  torus,  which  was  developed  and  coded  by  Wiesel.  In  addition,  the  new  theory 
is  a  complete  orbit  determination  method  which  retrieves  the  TLE  data  for  two  months  from 
www .  space-track .  com  servers,  propagates  each  TLE  for  one  orbit,  build  the  KAM  torus 
theory  file  using  some  of  the  SGP4  mean  elements,  and  fits  the  low  eccentricity  KAM  torus 
to  the  SGP4  and  TLE  predictions  using  least  squares  without  human  intervention.  Although 
converting  whole  TLE  catalog  to  the  new  method  is  a  trivial  task,  it  requires  at  least  30  days 
of  computation  time.  Therefore,  the  author  pseudo-randomly  chose  1500  objects  from  the 
NORAD  satellite  catalog  for  analysis.  This  chapter  can  simply  be  divided  into  four  dis¬ 
tinct  parts.  The  first  part,  Section  3.1,  discusses  the  data  gathering,  SGP4  propagation,  and 
contains  a  quick  overview  of  the  whole  process.  The  second  part  includes  all  sections  from 
Section  3.2  to  Section  3.6,  which  can  be  conveniently  named  as  the  KAM  theory  building 
part.  The  third  part,  starting  in  Section  3.7,  introduces  the  evolutionary  steps  to  the  desired 
KAM  torus  by  differential  correcting  the  initial  conditions  and  the  frequencies  of  the  newly 
built  torus  using  SGP4  and  TLE  data.  The  fourth  part  is  a  summary. 

3.1  Data  Gathering  and  Overview  of  the  Method 

This  effort  uses  the  TLE  obtained  from  the  website  www .  space-track .  org.  This  site 
allows  users  to  download  the  TLE  of  most  satellites.  The  previous  work  by  Wiesel  and  oth¬ 
ers  showed  that  there  are  some  properties  of  satellites  which  cause  difficulties  in  KAM  tori 
construction  process.  The  author  also  made  an  initial  test  in  order  to  define  an  optimal  ec- 
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centricity  threshold  level  because  this  effort  is  intended  for  low  eccentric  orbits.  Although 
the  author  defined  an  eccentricity  of  10  3  and  small  as  optimal  eccentricity  for  the  pro¬ 
gram  to  converge,  he  will  further  analyze  the  highest  limit  for  the  eccentricity  to  converge 
by  modifying  the  code.  It  is  also  known  that  there  are  no  nearly  circular  periodic  orbits  in 
the  close  vicinity  of  the  critical  inclination,  which  is  63.435°,  or  the  complementary  angle 
of  116.565°  [58].  Although  there  exist  periodic  orbits  for  polar  orbits,  the  close  vicinity 
of  polar  orbits  is  avoided  because  of  numerically  unstable  Legendre  polynomial  recursions 
for  the  gravitational  potential  calculations.  The  correction  is  left  for  future  studies.  An 
altitude  of  300  km  and  more  is  chosen  as  a  criterion  because  air  drag  shrinks  the  size  of  the 
torus  down  in  phase  space.  It  is  also  known  that  orbital  maneuvers  destroy  the  KAM  torus, 
but  they  are  not  listed  as  a  criterion  because  maneuvers  are  unpredictable.  Table  8  shows 
known  issues  in  the  previous  works  and  the  selection  criteria  of  test  objects  for  avoiding 
these  issues. 

Table  8.  Test  Case  Selection  Criteria 


Issues 

Test  Case  Selection  Criteria 

Air  Drag 

Objects  which  have  altitute  of  300 
kilometers  and  more 

Resonance 

Objects  with  periods  that  are  not 
nearly  an  integer  multiple  of  Earth’s 
rotational  period 

Critical  inclination  and 
Polar  Orbits 

Objects  which  have  inclinations  that 
aren’t  in  the  close  vicinity  of  the  criti¬ 
cal  and  polar  inclinations. 

(0°  <  i  <  60°  and  67°  <  i  <  87°  and 
93°  <  i  <  113°  and  119°  <  i  <  180°) 

There  are  7,938  space  objects  which  match  the  criteria  in  Table  8.  The  author  wrote  a 
C++  script  which  retrieves  7,938  TLEs  for  a  period  of  2  months  , which  is  between  Jan  1, 
2013  and  March  1,  2013,  from  www.space-track.org  server.  This  script  also  prepares 
essential  input  files  that  are  required  by  low  eccentricity  KAM  torus  construction  program, 
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TORIS  1 


Figure  16.  The  Orbit  of  Toris  1  Satellite  over  2  Months 

which  can  be  conveniently  named  as  theory  building  program.  The  solar  activity  is  rela¬ 
tively  low  between  the  selected  dates.  The  dates  for  solar  minimum  can  be  selected,  but 
the  author  intended  to  analyse  the  recently  launched  satellites  as  well.  Then,  1500  TLEs 
are  pseudo-randomly  selected,  and  propagated  for  an  orbit  for  a  period  of  2  months  by 
Vallado’s  revised  SGP4  code  [21].  Figure  16  shows  the  orbit  of  TORIS  1  satellite  over  two 
months,  which  is  simply  a  torus  shape.  TLEs  are  propagated  for  an  orbit  period  because  it 
is  known  that  the  position  errors  are  smallest  some  point  within  the  data  interval  and  grows 
with  respect  to  time  going  outwards  into  the  future,  or  back  into  the  past  [67].  It  is  also 
desirable  that  JSpOC  releases  the  TLE  with  an  epoch  time  within  data  interval.  JSpOC 
fits  the  orbits  for  3  to  4  days  for  LEO  satellites  and  a  couple  of  weeks  for  higher  altitude 
satellites.  Therefore,  there  is  a  progression  of  orbits,  and  that  progression  includes  more 
information  than  each  individual  orbit  has.  The  idea  behind  propagating  each  TLE  for  a  2 
months  orbit  is  to  obtain  more  accuracy  through  combining  smaller  chunks  of  less  accurate 
propagation  data.  Because  SGP4  and  TLE  prediction  data  are  considered  observational 
data,  data  from  the  calculation  of  dynamics  are  needed,  which  will  be  combined  to  form 
an  estimate.  The  low  eccentricity  KAM  torus  provides  the  dynamical  data.  The  low  ec¬ 
centricity  KAM  torus  construction  starts  with  the  calculation  of  the  periodic  orbit,  which 
is  basically  a  boundary  value  problem.  Once  the  periodic  orbit  is  built,  it  is  transformed 
to  Fourier  series  [23].  The  solution  for  the  periodic  orbit  in  the  zonal  potential  problem 
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yields  two  basis  frequencies,  which  are  the  Keplerian  frequency  and  the  nodal  regression 
rate.  The  other  missing  frequency,  which  is  the  advance  of  the  argument  of  perigee  rate, 
can  be  obtained  from  the  Floquet  solution  of  the  periodic  orbit.  The  Floquet  solution  also 
yields  two  time  linear  terms  due  to  adjacent  periodic  orbits.  The  solution  includes  only 
the  zonal  potential  so  far.  However,  perturbations,  which  are  air  drag,  and  sectoral  and 
tesseral  harmonics  perturbations,  are  added  to  the  Floquet  solution.  Therefore,  the  problem 
becomes  forced  linear  system  problem.  The  order  of  perturbations  suddenly  drops  from 
10  3  to  10  ,  which  satisfies  the  main  dictate  of  KAM  theorem  of  small  perturbations  bet¬ 

ter.  In  addition,  the  momenta  are  inertial  velocity  components  defined  in  a  rotating  frame 
of  reference  because  of  the  set  up  of  the  dynamics.  Therefore,  simply  substituting  veloc¬ 
ity  components  with  momenta  is  a  valid  approach,  see  Section  3.2.  Indeed,  defining  the 
forcing  function  as  a  function  of  Q\  and  Qi  in  a  reference  frame  that  rotates  with  Earth 
where  sectoral  and  tesseral  perturbations  are  stationary  leads  to  the  usual  set  up  of  the  fun¬ 
damental  equation  for  the  KAM  theorem,  see  Equation  86.  After  the  low  eccentricity  KAM 
torus  is  built,  it  is  fitted  to  the  observational  data  acquired  from  SGP4  and  TLE  prediction. 
Non-linear  least  square  algorithm  is  used  to  form  an  estimate.  The  modal  variables  and  the 
frequencies  of  the  periodic  orbit  are  updated  for  each  iteration,  see  Section  3.7.  Although, 
the  author  represents  the  process  piece  by  piece,  each  main  program  is  connected  to  another 
by  GNU/Linux  bash  files,  which  are  batch  files  for  MS-DOS  environment.  Therefore,  it  is 
a  simple  matter  to  switch  different  parameters,  such  as  air  drag,  and  sectoral  and  tesseral 
perturbations,  on  and  off.  Moreover,  human  intervention  is  not  required  during  the  pro¬ 
cess  due  to  the  bash  files  which  orchestrate  the  whole  orbit  determination  procedure.  This 
method  is  a  numerically  constructed  perturbation  theory  which  uses  periodic  orbits  under 
the  influence  of  only  the  zonal  potential  as  the  starting  point  [69]. 
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3.2  Satellite  Dynamics 


This  effort  uses  different  reference  frames,  such  as  rotating  and  inertial  frames  of  refer¬ 
ence.  The  inertial  reference  frame  is  the  ECI  frame,  and  the  rotating  frames  have  rotation 
rates  ,  (0  —  Cl  where  the  frame  of  reference  is  tied  to  the  regression  of  the  ascending  node, 
or  co  =  a)®  where  it  is  tied  to  Greenwich.  In  such  a  rotating  frame  of  reference  the  inertial 
velocity  vector  becomes  [69] : 


x  —  coy 

V  =  r=  <  y+CQx  > 


(96) 


where  (0  appears  in  the  components  of  the  inertial  velocity  because  of  transport  theorem. 
Then  the  satellite’s  kinetic  energy  becomes  [69]: 


T  =  \  {{x^COy)2  +  (y  +  (Ox)2  +  (z)2) 


(97) 


The  generalized  momenta  can  be  found  from  the  canonical  Hamilton  equations,  pi  = 
dT /dqi,  as  [69]: 


Px  —  x  —  (Oy 

py—y  +  cox  (98) 

Pz  =  z 
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where  the  generalized  momenta  are  the  inertial  velocity  components.  Then,  the  Hamilto¬ 
nian  can  be  calculated  from  H  =  ptq  —  T  +  V  as  [69] : 
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(99) 


where  /i  is  the  Earth’s  gravitational  parameter,  7?®  is  the  radius  of  Earth  at  the  equator,  and 
Cnm,  Snm  are  the  geopotential  coefficients.  This  effort  uses  NASA  EGM-96  gravity  model 
with  an  order  and  degree  of  20.  The  relationship  between  rectangular  coordinates  and  the 
radius,  r,  the  latitude,  8,  the  longitude,  A,  can  be  shown  as  [69]: 


r  =  \Zx2  +y2  +z2 


sin  8 
tan  A 


z 

\Jx2+y2 

y_ 

X 


(100) 


where  the  rectangular  coordinates  are  in  a  reference  frame  tied  to  the  Greenwich,  which 
is  the  Greenwich  referenced  ECR  frame,  see  Section  2.4.  The  expansion  of  the  geopo¬ 
tential  in  the  zonal  harmonics,  which  is  m  =  0,  is  used  for  constructing  the  periodic  orbit. 
The  sectoral  and  tesseral  harmonics  are  added  as  perturbations  to  the  periodic  orbit.  The 
Hamilton’s  equations  can  be  represented  as  [69]: 


X 


(101) 


where  X  is  the  physical  state  vector,  XT  —  {.r.y. z,  px,  py.  /?-  },  and  Z  is  a  2 n  x  In  matrix, 
where  n  is  the  number  of  coordinates,  that  conveniently  arranges  the  Hamilton’s  equations 
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in  vector  form  as  [69]: 
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(102) 


The  linearization  of  the  Hamilton’s  equations  is  needed  because  this  effort  is  an  estimation 
problem.  The  equations  of  variation  can  be  represented  as  [69]: 


d2H 

i  =  Z3X2=Ax’ 


(103) 


where  x  is  represented  as  small  changes  to  the  physical  state  vector  X. 

Air  drag  perturbations  are  also  added  to  the  the  combination  periodic  orbit  and  the 
Floquet  solution.  The  momenta  are  inertial  velocity  components.  Therefore,  air  drag  ac¬ 
celeration  can  be  represented  as  [69]: 


_  '  CdApo  P  /  2  |  „2  |  „2n 

admg  ~  ~2  m  V  Px +py  +  Pz  P 
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Po 
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(104) 


where  B*  is  another  measure  for  the  susceptibility  of  a  satellite  to  air  drag  ,  m  is  the  satellite 
mass,  and  p  is  the  atmospheric  density  at  perigee.  [60].  For  further  information,  see 
Equation  84. 
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Symmetric  About  Polar  Axis 


Figure  17.  The  Earth’s  Zonal  Gravity  Potential 


3.3  Periodic  Orbits 

The  periodic  orbit  and  Floquet  theory  are  used  as  starting  point  for  calculating  the 
perturbations  because  the  previous  KAM  torus  construction  methods  failed  for  low  eccen¬ 
tricity  orbits  and  the  periodic  orbit  resides  at  the  core  of  the  torus.  The  modified  Laskar 
algorithm,  developed  by  Wiesel,  is  one  of  the  KAM  torus  construction  methods.  It  yields 
three  basis  frequencies  for  orbits  by  Fourier  series  spectral  analysis.  The  low  eccentricities 
are  hard  to  be  detected  for  methods  like  Laskar  frequency  analysis  method,  which  requires 
high  eccentricities  that  cause  more  dramatic  change  in  the  frequency  spectrum  of  the  orbit. 
Although  this  fact  is  a  disadvantage  for  the  previous  KAM  torus  construction  methods,  it  is 
favorable  for  this  effort  because  the  structure  lies  at  the  core  of  the  torus  is  a  periodic  orbit. 

In  1960’s,  the  periodic  orbits  were  studied  extensively  because  they  are  more  realistic 
than  2BP,  which  describes  every  orbit  periodic.  The  Newtonian  point  mass  plus  the  zonal 
harmonics  terms  lead  to  the  periodic  orbit  because  the  zonal  potential  about  polar  axis  is 
symmetric.  Figure  17  represents  the  symmetry  of  the  Earth’s  zonal  potential  about  polar 
axis.  Nearly  circular  orbits  are  periodic  orbits  in  the  Earth’s  zonal  gravity  potential  for 
every  inclination,  except  in  the  close  vicinity  of  the  critical  inclination  [58].  The  periodic 
orbits  are  closely  related  to  the  frozen  orbits,  which  have  nearly  stationary  eccentricity 
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n 

Figure  18.  The  Regression  of  Orbit  Plane  in  the  Zonal  Gravity  Potential 

and  argument  of  perigee.  The  repeated  ground  tracks  are  of  great  importance  for  mission 
planning.  The  periodic  orbits  defined  in  a  reference  frame  tied  to  the  rotating  Earth  have 
repeated  ground  tracks.  Therefore,  the  periodic  orbits  have  been  studied  as  frozen  orbits 
as  well  [44,  43].  It  has  been  shown  that  nearly  circular  orbits  are  unstable  at  the  criti¬ 
cal  inclination,  and  the  periodic  orbits  in  the  vicinity  of  the  critical  inclination  have  high 
eccentricities  [17].  The  numerical  proof  of  the  existence  of  periodic  orbits  with  high  ec¬ 
centricities  has  been  done  by  Broucke  [8].  Moreover,  Wiesel  has  developed  a  new  relative 
motion  solution  for  satellites  using  the  low  eccentricity  periodic  orbits,  but  their  potential 
use  for  orbit  determination  hasn’t  been  explored  until  he  uses  the  low  eccentricity  periodic 
orbits  to  construct  the  low  eccentricity  KAM  torus  [69] . 

In  the  zonal  problem,  there  are  periodic  orbits.  However,  there  are  a  few  difficulties  in 
building  the  periodic  orbit.  It  is  known  that  there  are  two  non-zero  frequencies,  which  are 
the  Keplerian  frequency  and  the  regression  of  the  node.  These  frequencies  aren’t  known. 
The  periodic  orbit  precesses  at  a  certain  rate  which  can  only  be  determined  when  it  fin¬ 
ishes  its  motion  for  one  period,  and  the  orbits  aren’t  periodic  in  inertial  space.  Therefore, 
defining  the  periodic  orbit  in  a  reference  frame  that  rotates  with  the  node  is  required.  Fig- 
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ure  18  represents  the  regression  of  orbit  plane  under  the  influence  of  zonal  potential.  The 
periodicity  condition  needs  to  be  set  up  in  such  a  way  that  the  nodal  regression  rate  of  the 
orbit  plane  is  not  needed.  Therefore,  the  initial  conditions  are  defined  with  position  vector, 
r  =  {xo.0.0},  and  momentum  vector,  p  =  {x0,v0cos  i0,v0sin  ia}.  The  orbital  period,  T,  and 
inclination,  i,  define  the  orbit,  while  the  unknowns  are  ZT  =  {x0,x0,  va}.  Because  Earth 
doesn’t  have  North/South  symmetry,  xa  is  required  as  a  component  of  the  momentum  vec¬ 
tor,  p.  The  orbit  becomes  periodic  if  the  3  conditions  after  one  revolution,  t,  are  set  up  as 
[69]: 

z(t) 

G=<  r(T)-p(t )-x0x0  7  =  0,  (105) 

r(T)-r(r  )-x20 


where  the  conditions  dictate  the  satellite  to  cross  the  equatorial  plane,  z(t)  =0,  to  have 
the  same  radial  velocity  as  at  the  beginning,  r(r)  -p(r)  —x0x0  —  0,  and  to  be  at  the  same 
distance  from  Earth  as  at  the  beginning.  The  actual  orbit  precession  isn’t  known,  but  can 
be  calculates  as  [69]: 


1  _i  f(t)  ■  r(0)  1  _!*«) 

-COS  ,  7  —  /  ,  =  -COS  — — 

r(T)|jr(0)|  T  *(0) 


(106) 


A  linearization  of  Equation  105  using  Taylor’s  series  expansion  is  required  to  obtain 
the  periodic  orbit  [69] : 

G  =  G(S)  +  —  5S  =  0  (107) 

do. 

The  expansion  using  the  state  transition  matrix  T>(t,0)  ,  which  is  calculated  from  parallel 
numerical  integration  with  the  trajectory,  yields  the  derivative  matrix  of  G  [69]: 
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where  T>(t,0)  can  be  calculated  from  <9X(t)/<9X(0).  The  partial  derivatives  in  the  Equa¬ 
tion  107  can  be  computed  as  [69]: 
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i  y  z  x  y  z  > 

2x  2 y  2z  0  0  0 

> 

0  0  0  0  0  0 

—xQ  0  0  —x0  0  0  > 

— 2x0  0  0  0  0  0 


dx(0) 

dZ 


10  0  0 
0  0  0  1 
0  0  0  0 


0  0 

0  0 

cos  i0  sin  iQ 


(109) 


Then,  the  initial  condition  parameter  vector,  E,  can  be  corrected  by  the  correction,  53, 
which  can  be  calculated  by  Equation  107  [69]. 

After  the  periodic  orbit  is  built,  it  is  represented  as  a  Fourier  series  by  harmonic  analysis 
[23].  Figure  19  shows  the  representation  of  the  periodic  orbit  by  two  angle  variables.  In 
Figure  19,  Q\  is  ’’mean  argument  of  latitude  analogue”,  and  Qi  is  regression  of  the  right 
ascension  of  ascending  node.  Q\  and  Qi  are  the  physical  coordinates  of  the  KAM  torus 
[69]: 

Ql=M  +  a,  (110) 


Q2  =  n-dg,  (111) 

where  M  is  mean  anomaly,  Q  is  right  ascension  of  ascending  node,  and  0„  is  angle  between 
the  vernal  equinox  and  the  prime  meridian.  The  periodic  orbit  is  a  function  of  only  Q\  in  a 
frame  of  reference,  which  is  the  ECNF  frame,  see  Section  2.4,  that  rotates  with  the  ascend- 
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Figure  19.  The  Representation  of  the  Periodic  Orbit  by  Angle  Variables 


ing  node.  The  nodal  orbital  frequency,  6)\,  and  the  inertial  nodal  regression  rate,  Oh,  can 
be  obtained  from  the  Fourier  series.  These  Fourier  series  store  only  the  coordinates,  ,  77, 
and  £,  because  the  momenta  can  be  obtained  by  the  kinematic  relations  in  Equation  98. 
Moreover,  it  is  a  simple  matter  to  generalize  from  ECEF  frame  to  ECI  or  ECEF  frames 
because  the  zonal  potential  is  symmetric  about  polar  axis.  The  inertial  nodal  regression 
rate,  Oh,  needs  to  be  updated  to  define  the  regression  rate  with  respect  to  the  Greenwich. 
Equation  112  represents  that  the  periodic  orbit  defined  in  ECNF  frame,  where  the  periodic 
orbit,  X,  is  a  function  of  only  <2t,  can  be  transformed  to  ECI  or  ECEF  frame  with  a  double 
rotation  matrix  about  the  z  axis  by  the  nodal  angle  Q2  [69]: 


Rzi-Qi)  0 
0  Rz{-Qi) 


(112) 


Equation  112  can  also  be  used  when  the  desired  reference  frame  is  ECI  instead  of  ECEF 
with  small  modifications.  Both  Qi  in  the  equation  are  substituted  by  the  usual  initial  frame 
node,  Q.,  and  Oh  is  reduced  by  the  Earth’s  rotation  to  transform  the  periodic  orbit,  X,  from 
ECNF  frame  to  ECI  frame.  Although  the  two  frequencies  have  been  obtained  from  the 
periodic  orbit,  the  third  frequency  is  required  for  the  solution. 


67 


3.4  Floquet  Solution 


Periodic  orbits  plus  the  Floquet  solution  provides  a  method  that  numerically  solves  dy¬ 
namical  problems  which  don’t  have  closed  form  solution.  The  linearization  of  the  periodic 
orbit  describes  the  behaviour  of  the  system  near  the  periodic  orbit.  Linearizing  about  the 
periodic  orbit  results  in  a  linearization  problem  [64]: 

x  =  A{t)x ,  (113) 

where  x  is  the  first  order  displacements  from  the  state  X.  Equation  113  is  a  time  periodic 
linear  differential  equations.  The  solution  to  time  periodic  linear  differential  equations 
yields  stability  information  of  periodic  orbits.  The  solution  to  the  Equation  113  can  be 
obtained  as  [64]: 

®{t,t0)  =  E(t)exp(J(t -to))E~\to),  (114) 

where  £  is  a  periodic  matrix,  /  is  a  matrix  of  system  frequencies  in  the  Jordan  normal  form. 
The  system  frequencies  are  also  called  Poincare  exponents. 

The  Floquet  solution  can  be  obtained  by  calculating  the  Jordan  normal  form,  J,  and  the 
periodic  matrix,  E,  over  one  period.  Over  one  period,  Equation  114  becomes  [64]: 

<J>(t,0)  =E(t)exp(Jt)E~l(0),  (115) 

where  E(t)  =  E( 0)  because  E  is  periodic.  Figure  20  represents  the  visual  depiction  of 
E  periodic  matrix.  Geometrically,  defining  E{t)  —  E( 0)  puts  a  sets  of  coordinates  on  the 
periodic  orbit,  and  they  follow  the  trajectory  by  rotating  around  themselves.  At  the  end 
these  sets  of  coordinates  smoothly  join  themselves.  <4>  is  the  state  transition  matrix  over  a 
period.  Indeed,  because  the  system  is  periodic,  we  have  <t>  for  all  time.  The  <t>  matrix  is 
called  monodromy  matrix.  The  eigenvalues  and  eigenvectors  of  the  <t>  matrix  can  be  found 
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Figure  20.  E  Periodic  Matrix 


by  rearranging  Equation  114  [64]: 


£(0)_1<I>(r,0)E(0)  =eJt, 


(116) 


where  the  exponential  of  the  Jordan  form  ,J,  is  the  eigenvalue  matrix  of  <t>.  E  periodic 


matrix  yields  the  eigenvalues  of  the  <t>  matrix.  If  A/  are  eigenvalues  of  the  monodromy 
matrix,  Equation  116  becomes  [64]: 


(117) 


where  ft),  are  the  system  frequencies,  which  are  also  called  Poincare  exponents.  The  inter¬ 
pretation  of  the  Poincare  exponents  is  similar  to  the  eigenvalues  of  a  constant  coefficient 
system.  For  Hamiltonian  systems  without  dissipation,  the  frequencies  must  be  a  pair  of 
positive  /  negative  real  numbers  or  a  pair  of  a  pair  of  imaginary  numbers.  As  an  exception, 
when  A,  =  0,  there  is  a  repeated  value  of  ft)/  =  0.  Repeated  roots  in  linear  systems  lead  to 
off-diagonal  time  terms.  These  pair  of  zeros  are  called  a  degenerate  mode.  For  this  effort, 
the  Hamiltonian  has  two  integrals  of  the  motion,  which  are  the  Hamiltonian  itself ,  and  the 
z  component  of  the  angular  momentum.  Because  each  of  these  pair  of  zeros  correspond 
to  one  integral  of  motion,  this  effort  has  two  degenerate  modes.  Each  degenerate  mode 
involves  a  static  displacement  and  a  linear  drift  is  associate  momentum  is  not  zero.  These 
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static  displacements  can  occur  either  along  the  orbit  or  around  the  polar  axis.  These  two 
degenerate  modes  correspond  to  two  frequencies.  The  third  frequency  is  the  free  eccen¬ 
tricity  of  the  orbit,  which  is  a  pair  of  imaginary  frequencies  by  the  eigenvalue  calculation, 
because  the  imaginary  part  of  ft),  is  the  oscillatory  frequency  of  the  mode  i.  For  this  effort, 
the  oscillatory  frequency  is  ±©3  with  algebraic  sign  that  depends  on  the  inclination.  As  a 
result,  the  calculation  of  eigenvalues  yields  the  complete  linear  stability  information  of  the 
periodic  orbit  [64,  69]. 

The  calculation  of  eigenvector  yields  E(t)  over  one  period,  which  can  be  used  to  trans¬ 
form  from  physical  variables  to  model  variables.  The  eigenvalue  problem  for  equation  is 
highly  singular,  and  a  special  treatment  is  required.  The  differentiation  of  the  state  transi¬ 
tion  matrix  becomes  [63]: 


<b{t1t0)=E(t)eJit-to)E~l(t0)+E(t)JeJ{t-to)E-l(t0)  (118) 


By  substituting  into  Equation  119  : 

<P(t,t0)=A(t)<P(t,t0)  (119) 

A  differential  equation  for  F(t)  becomes: 

E(t)eJ(t-to)E-\t0)+E(t)JeJ(t-t°)E-\t0)  =  A(t)E(t)exp(J(t  - t0))E~l  (t0) 

E(t)  +E(t)J  =A(t)E(t)  (120) 

E  —AE  —  EJ 

Since  E(0)  is  obtained  from  the  periodic  orbit,  solving  Equation  120  from  /  =  0  to  /  yields 
F(t)  over  one  period.  E(t  )  periodic  matrix  is  reduced  to  Fourier  series  [23].  For  this  effort, 
the  periodic  matrix  is  a  function  of  Q\,  E{Q\). 

The  transformation  from  physical  coordinates  to  modal  variables  is  possible  because 


70 


F(t)  periodic  matrix  is  calculated.  Assuming  y  is  the  modal  variable  [63]: 


y  =  E~l(t)x 


(121) 


By  substituting  into  Equation  122  : 


x(t)  =  (j>(t,t0)x(t0) 


(122) 


The  Floquet  solution  represented  in  modal  variables  becomes  : 

Table  9.  Analogous  of  Modal  Variables 


Modal 

Variables 

Analogous  of  Modal  Variables 

yi 

a  displacement  along  the  argument  of  latitude  di¬ 
rection,  <2i- 

yi 

a  change  in  the  momentum  which  is  analogous  to 
the  Delaunay  momentum,  L. 

T3 

an  offset  in  the  nodal  variable,  Qi. 

T4 

a  change  in  the  z  component  of  orbital  angular  mo¬ 
mentum. 

ecos  (0 

ye 

esin  co 

x(t)=E(t)eJ^E~\t0)x(t0) 

E~l{t)x{t)  =  eJ^~to^E~l{t0)x{t0)  (123) 

y(t)  =  eJ[t-to)y{to) 


Equation  123  is  a  solution  to  the  linear  system  to  the  form  [69]: 


y  =  Jy 


(124) 
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The  Floquet  solution  can  be  added  to  the  periodic  orbit  solution  as  : 


*ECR  =  Rz(-Qi)  {XPO(Qi)  +  E(Qi)exp(J(t  -  t0))y(t0)} ,  (125) 


where  double  z  rotation  is  represented  as  R-.  the  nodal  frame  periodic  orbit  series  are  shown 
as  (. X)po .  The  modal  variables  ,  y\  and  yi,  are  local  versions  of  the  global  angles  ,  Q\  and 
Q2 ,  and  they  are  scaled  to  suit  this  purpose.  The  Jordan  form  matrix  becomes  : 


exp(J(t  -t0))  =  < 


1  t  —  t0  0  0 

0  10  0 

0  0  1  t  —  ta 

0  0  0  1 

0  0  0  0 

0  0  0  0 


0  0 

0  0 

0  0 

> 

0  0 

cos  <23  sin  Qi 
—sin  <23  cos  Q3 


(126) 


The  modal  variables  are  analogous  to  6  orbital  elements  that  can  be  used  to  represent 
nearby  motion.  Table  9  represents  an  analogy  between  modal  variables  and  miscellaneous 
orbital  elements.  The  global  variables,  <2i  and  Qi,  should  be  freely  traded  with  the  modal 
variables,  y\  and  y$,  because  this  will  keep  them  from  growing  unboundedly  [69]. 


3.5  First-Order  Perturbations 

The  periodic  orbit  and  its  Floquet  solution  has  been  obtained  from  Equation  125.  This 
solution  has  the  Newtonian  term  and  all  of  the  zonal  gravity  harmonics.  The  perturbations 
start  approximately  on  the  order  of  10  5  instead  of  1 0  3 ,  which  is  advantageous  because 
errors  originated  from  the  calculation  of  perturbations  are  less  influential.  However,  the 
method  for  calculating  first-order  perturbations  can  be  extended  to  higher  order  perturba¬ 
tions.  This  effort  includes  only  sectoral  and  tesseral  harmonics  and  air  drag.  However,  the 
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exact  same  method  applies  to  lunar  and  solar  mass  perturbations  by  including  at  least  one 
extra  angle  in  the  Fourier  series  which  describes  the  motion  of  the  Moon  or  the  Sun  with 
respect  to  Earth  in  a  frame  that  rotates  with  it.  The  perturbing  force,  Xpert  ,  can  be  added  to 
Equation  113  as  [69]: 

X  =  Ax  +  Xpert  ( 1 27) 


The  perturbing  acceleration  can  be  expanded  about  the  periodic  orbit  because  this  is  a 
forced  linear  system  [69] : 


X 


X 

=  X 


pert  ^  ^pert\Xpo+  ^ 

dx 


dx 


pert 


X+  ... 


X 


PO 


pert|Xpo+  dX 


pert 


(128) 


Ey  + ... 


Xpo 


The  first-order  perturbations  can  be  added  to  the  Floquet  solution  by  Equation  124  as 
[69]: 

(  \ 


1  0  0  0  0  0 
0  1  0  0  0  0 


3>  = 


0  0  10 
0  0  0  1 


0 

0 


0 

0 


}y  +  E-lX 


perty 


0  0  0  0  0  003 

0  0  0  0  — ft)3  0 


(129) 


where  for  sectoral  and  tesseral  perturbations,  the  trailing  term,  E~lXpert,  is  a  function  of 
the  global  variables,  Q\  and  (?2-  The  forcing  term  can  be  transformed  to  a  double  Fourier 
series  in  these  angles.  This  applies  to  air  drag  as  well.  However,  the  air  drag  forcing  term 
is  a  function  of  only  one  global  variable,  Q\. 

For  the  first-order  perturbations,  a  perturbing  acceleration  transformed  to  the  modal 
variables  is  decoupled.  Two  different  terms,  which  are  the  constant  and  periodic  Fourier 
series  terms,  are  examined  for  the  perturbing  acceleration  as  the  response  to  the  oscilla- 
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tory  and  the  two  degenerate  modes.  Each  factor  is  evaluated  separately.  A  constant  solu¬ 
tion  is  obtained  from  the  constant  terms  in  the  Fourier  series  for  the  oscillatory  mode  by 
Equation  129  as  [69]: 


/ .  \ 

f  A 

_ 

0 

{y6) 

(°) 

0  (O3 

-(Os  0 


(» 


M 


w  \C6J 


Then,  the  solution  to  the  Equation  130  becomes  [69]: 


yA 

(  C6/ (O3 

* 

^  -c5/(03 

(130) 


(131) 


These  terms  can’t  be  ignored  because  (O3  is  a  small  number  for  the  most  part.  The  modal 
variables,  y\  and  >’2,  don’t  become  infinite  near  the  periodic  orbit,  which  has  zero  eccen¬ 
tricity  for  this  effort,  because  they  are  analogous  to  the  quantities,  ecos  (O  and  esin  00,  in 
the  2BP.  The  constant  terms  in  the  Fourier  series  for  one  of  the  degenerate  modes  doesn’t 
yield  a  constant  solution  [69] : 


/  .  \ 

yi 

w 


0  1 
0  0 


(n 


\ 


\ 


A2)  \C2J 


Then,  the  solution  becomes: 


l  \ 

/ 

\yi 

_ 

\y2J 

V 

(ci  -  y2{t0))(t-to)  +  \c2{t 2 
C2t-yi{t0) 


) 


(132) 


(133) 


where  the  solution  for  3)2  =  C2  is  written  at  the  start  time  because  the  other  secular  terms, 
which  appear  due  to  the  periodic  perturbations  in  y2,  can  be  absorbed.  The  solution  for  y\ 
has  been  obtained  by  substituting  V2  into  y\.  Equation  133  represents  that  y  1  and  >’2  could 
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have  quadratic  terms,  and  y2  and  34  could  have  a  linear  rate  of  change.  When  there  is  no 
dissipative  force  in  the  system,  C2  should  be  zero,  and  a  quadratic  term  isn’t  expected  for  in¬ 
track  and  nodal  perturbation.  For  air  drag  perturbation,  quadratic  behavior  is  cxpcctcd.  vi 
and  y2  can  be  observed  by  their  associated  global  variables  because  the  modal  variables 
aren’t  scaled.  The  exact  same  approach  is  taken  for  periodic  Fourier  terms. 

The  periodic  terms  in  the  Fourier  series  for  a  degenerate  mode  yields  the  forced  linear 
system  as  [69]: 


l  •  \ 

y\ 

w 


\  0  0  J  U  y 


( 

+ 

V 


Cl  cos{n\Qx  +n2Q2)  +si  sin(n\Q\  +  n2Q2) 
c2  cos(niQi+n2Q2)  +s2  sin{niQ\  +  n2Q2) 


\ 

) 


(134) 


where  the  setup  is  for  sectoral  and  tesseral  harmonics,  but  if  a  perturbation  requires  it, 
an  additional  angle  can  be  introduced.  In  the  equation,  n\  and  n2  are  integers,  and  wp  = 
n\(D\  +n2ah.  The  solution  for  y2  can  be  calculated  [69]: 


yi(t)  =  ( —  sin(niQi  +n2Q2) - cos(niQ\  +n2Q2 )) 

(0p  (Op 

Then,  substituting  into  the  yi  yields  : 


(135) 


yiO) 


^1 

sin(niQi  +  n2Q2)  -  ( — 
(Op 


)  cos{rnQi+n2Q2) 


C2  S  2 

(—  sin{n\Q\  +n2Q2 ) - -  cos(niQi  +n2Q2)) 

COp  (Op 


t 


to 


(136) 


where  the  constant  terms  in  the  V2  solution  lead  to  time  dependent  terms  in  the  solution 
for  y  1 .  However,  the  secular  terms  introduced  by  the  constant  terms  in  the  V2  solution  has 
already  been  absorbed  into  the  solution  provided  by  Equation  133. 
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The  remaining  forced  linear  system  is  due  to  the  periodic  terms  in  the  oscillatory  mode. 


Equation  137  represents  the  forced  linear  system  [69]: 


(  . 

ys 

\ ye 


\ 

/ 


J  0  0)3 

j  -®3  0 


(  \ 

ys 

\y6J 


c5  cos{rnQi+n2Q2)  +S5  sin(mQi  +  n2Q2) 
c6  cos(niQi  +n2Q2)  +  s6  sin{n\Q\  +n2Q2) 


\ 

/ 


Let  both  a  forced  solution  and  the  perturbation  term  has  the  same  form: 


(  \ 

ys 

W 


a5  cos(niQi+n2Q2)  +  p5  sin{n\Qi  +  n2Q2) 

^  a6cos(niQi+n2Q2)+P6sin(niQi+n2Q2)  ^ 


(137) 


(138) 


After  inserting  assumed  solution  and  simplifying  it,  the  solution  becomes  a  set  of  linear 
equations  for  the  forcing  coefficients  as: 


as 

a6 


C6(03+S5(0p 

®32-® i 

-C5(03+S6(0p 
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S(,(d3—C5(Op 
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—S50>3  ~  C6(Op 
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(Of 


(139) 


Once  the  perturbing  Fourier  series  is  obtained,  it  is  converted  to  the  perturbation  solution, 
and  stored.  This  method  can  be  applied  to  any  other  perturbations  with  simple  modifica¬ 
tions,  and  the  code  written  for  this  method  requires  only  the  new  data  file  and  truncation 
level  associated  with  the  perturbation.  It  should  be  noted  that  Equation  139  includes  the 
expected  small  divisors  problem,  which  will  lead  to  problems  near  the  resonance. 

Equation  125  can  be  modified  to  include  the  first-order  perturbations  as[69]: 


XeCR  =Rz{-Ql)  {Xpo(2i)  +E(Q\)  [exp(J (t  —  t0))y free{t0)  +  }’ forced]  }  ,  (140) 
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where  the  modal  variables  from  the  Floquet  solution  is  represented  by  }'free,  and  the  per¬ 
turbed  solution  is  denoted  by  yforced ■  The  perturbed  solution,  y forced ■>  can  be  a  function  of 
different  angles,  depending  on  the  perturbation. 


3.6  Second-Order  Eccentricity  Perturbations 


The  modal  variables,  yi  and  y$,  are  necessary  for  fitting  orbit  and  for  small  perturba¬ 
tions.  The  global  angles  Q\  and  Qi  absorb  the  large  changes.  Moreover,  secular  terms  may 
be  absorbed  in  their  frequencies,  ft) i  and  ft>2.  However,  the  expansion  of  the  torus  along,  >’5 
and  yg,  which  describe  the  eccentricity  and  argument  of  perigee  states,  is  necessary  because 
not  all  orbits  have  small  enough  eccentricities.  The  equation  for  the  expansion  of  the  two 
model  variables,  y$  and  y^,  when  setting  the  other  model  variables  approximately  zero  in 
physical  variables  becomes  [69] : 


Xi  -  Zi 


d2H 


,adXadXp 


1 

x[3  +  ~  Zj 


d3H 


ia 


Xpo 

1 


dXadXpdXr 


Xfi  Xy  +  ... 


Xpo 


=  Aia(Q\)xa  +  xBjap  (Q\)xaxp 


(141) 


where  Z  is  the  symplectic  matrix,  i  =  (1...6)  for  this  specific  case,  the  linear,  time  varying 
differential  equations  are  the  equations  of  variation.  Representing  the  equation  in  the  modal 
variables,  Equation  141  becomes: 


1 

yi  =  Jiayoc  +  xEia  BapyEp8Ey£ysy£  +  ... 

f  (142) 

=  Jiaya  +  2B'iapy<xyp  +  ••• 

where  B\-k,  the  periodic  orbit,  and  the  E  matrix  are  periodic  functions  of  Q\.  The  quadratic 
matrix,  B'ijk  with  y$  and  y (,  portion  in  the  final  two  indices  yields  a  periodic  six  by  two  by 
two  tensor  [69]. 

Assuming  all  model  variables  are  zero  except  the  eccentricity  mode,  the  linear  matrix 
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becomes: 


Then,  the  linear  solution  can  easily  be  obtained  as: 


(143) 


y(0  =&(t,to)y(t0)  = 


cosOh,  (t  —  ta)  sindh,(t  — 10) 
—sinGh,[t  —  t0)  cosdh,{t  —  t0) 


(144) 


If  the  quadratic  terms  in  the  Equation  142  is  small  enough,  the  first  order  solution,  y(t)  = 
<t>{t,t0)y(t0),  can  be  used  for  the  equation  as: 

ji  =  Jiaya  +  ^'iaR^ay^  figyyityy  s(to) 

f  (145) 

=  Jiaya  +  ^’iapyaitobp  (*o) 

where  5”  is  a  function  of  both  Q\  and  Q 3  =  C03(t  —  to).  The  solution  to  Equation  145  is 
assumed  to  be  generalized  as  [69]: 


yt(t)  =  ®t?ya(t0)  +  ^ %ya{to)yp(t0 )  + ...  ,  (146) 

where  4>^l(t0,r0)  =  7,as  it  is  in  the  classical  linear  system  theory,  and  &2\t0,t0)  —  0  be¬ 
cause  the  initial  conditions,  yi(t0),  need  to  be  preserved.  Therefore,  the  time  derivative 
of  Equation  146  can  easily  be  calculated  because  yi{t0)  are  constant.  Inserting  the  time 
derivative  in  the  form  of  the  modal  differential  Equation  145  yields  [69]: 

yi  =  ^ya{to)  +  \^%ya{to)y^t0)  + ... 

=  Jip^plyM  +  \J^{yo.pya{to)yp  {to)  (147) 

2"®  iapya{to)yp{to) 
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where  the  arbitrary  indices  are  for  a  uniformity  requirement  in  the  initial  conditions.  The 
first  order  state  transition  matrix  differential  equation  can  be  obtained  as  [69]: 

(148) 

Then,  the  new  result  in  the  second  order  becomes: 

(149) 

(2) 

where  it  is  already  known  that  (t0,t0)  =  0.  The  solution  for  Equation  149  has  been 
provided  for  the  constant  terms  in  the  Fourier  series  of  5”  by  Equation  131  and  by  Equa¬ 
tion  139  for  the  periodic  terms.  If  indices  j  and  k  are  constant,  the  solution  to  the  system 
will  be  same  as  the  first  order  case.  Therefore,  <t>!2)  will  be  a  periodic  function  of  the  global 
variables,  Q\  and  Q^.  The  extension  to  higher  order  of  the  eccentricity /argument  of  perigee 
is  believed  to  better  characterize  the  short  periodic  oscillations. 

3.7  Least  Squares  Fit  to  SGP4  and  TLE  Data 

For  this  method,  the  point  on  the  periodic  orbit  that  corresponds  to  the  epoch  time  of 
SGP4  and  TLE  data  should  be  calculated  in  terms  of  Q\  and  Qi.  From  the  perspective  of 
estimation  theory,  a  reference  trajectory  should  be  corrected  to  yield  an  estimate  of  the  ac¬ 
tual  trajectory.  For  this  effort,  the  low  eccentricity  KAM  torus  predictions  are  considered  as 
the  reference  trajectory,  and  the  SGP4  and  TLE  predictions  are  considered  as  the  observa¬ 
tional  data.  Wiesel  and  Frey  showed  that  the  KAM  torus  basis  frequencies  can  be  obtained 
from  SGP4  and  TLE  data.  The  desired  periodic  orbit  is  the  one  with  zero  y 2  and  >’4  model 
variables.  The  least  squares  routine  for  this  effort  fits  the  low  KAM  torus  to  the  SGP4  and 
TLE  data  by  updating  the  global  variables,  Q\  and  Q2,  the  frequencies,  (0\  and  ah,  and  B*. 
The  convergence  criterion  is  met  when  the  estimate  correction  vectors  of  modal  variables, 
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Error  Ellipsoid 

Figure  21.  The  Converged  and  Unconverged  State  Correction  Vectors 

8yi,  are  smaller  than  1%  of  the  sigma  values,  O';.  Convergence  for  least  squares  method 
depends  on  the  correction  vectors  and  the  error  ellipsoid.  If  the  correction  vectors  lie  well 
within  the  1  o  error  ellipsoid  [67].  Figure  21  represents  the  converged  and  unconverged 
state  correction  vectors. 

This  effort  uses  nonlinear  least  squares.  This  method  combines  the  SGP4  and  TLE 
data,  which  has  been  propagated  for  an  orbit  over  two  months,  and  the  low  eccentricity 
KAM  torus  to  form  an  estimate.  For  each  observation,  6,  a  series  of  calculations  are  made 
for  nonlinear  least  squares  method.  First,  the  state  transition  matrix  is  obtained  , 
by  propagating  the  state  vector  to  the  observation  time,  f,-.  Second,  the  residual  vector, 
r,  =  Zi  —  G(x),  is  calculated.  The  SGP4  and  TLE  data  are  represented  as  zu  and  the  low 
KAM  torus  data  are  represented  as  G(jc).  Third,  the  observation  matrix,  T,  =  //,<$>,  the 
covariance  matrix,  P§x,  and  the  state  correction  vector  at  epoch  are  calculated  as  [67]: 

T,  =  H& 

Ps,  =  CLtTqt't,)-'  (150) 

i 

Sx(to )  =  Pdx^T^Q^n 

i 

where  Qi  is  a  matrix  that  has  sigma  errors  associated  with  the  measurements,  is  obser¬ 
vation  function  linearization  matrix.  Then,  the  reference  trajectory  is  corrected  as: 

Xref+l  do)  —  Xrefdo)  T  8x(to)  (151) 
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Table  10.  Three  Shortcomings  of  the  Least  Squares  Method  [11] 


1 

Although  every  observation  has  different  accuracy  for  the 
observations,  each  observation  is  weighted  equally,  which  is 
a  problem  when  there  is  observation  with  poor  accuracy  be¬ 
cause  the  least  squares  method  will  be  biased  towards  them. 

2 

The  errors  can  be  correlated,  not  independent  as  least 
squares  method  considers  them  to  be. 

3 

The  method  soesn’t  consider  that  the  errors  are  samples 
from  a  random  process,  and  it  doesn’t  utilize  any  statisti¬ 
cal  information 

A 


NULL  SPACE  OF  MATRIX  A  COLUMN  SPACE  OF  MATRIX  A 

Figure  22.  The  Visual  Depiction  of  the  Least  Squares  Method 


Finally,  the  process  ends  when  a  desired  limit  for  the  convergence  is  met.  However,  there 
are  short  comings  of  the  least  squares  solution,  which  are  listed  in  Table  10.  In  addition, 
Figure  22  shows  the  visual  depiction  of  least  squares  method.  The  least  squares  method 
yields  the  the  projection  of  the  solution  in  the  column  space,  where  the  solution  is  supposed 
to  exist.  Therefore,  in  a  sense,  the  least  square  method  minimizes  ( Ax  —  b )2,  and  yields  the 
closest  approximate  of  the  solution. 


3.8  Summary 

The  geometry  of  the  core  of  torus  in  the  full  geopotential  appears  to  be  a  multiply- 
periodic  two  dimensional  ribbon  structure  as  a  function  of  the  argument  of  latitude  angle, 
<2i,  and  the  nodal  angle,  Q2,  whereas  it  is  a  periodic  orbit  in  the  zonal  potential.  The  orbital 
eccentricity  is  a  mode  of  oscillation,  which  has  been  expanded  to  the  second  order.  For 
the  low  eccentricity  KAM  torus  theory,  each  satellite  has  a  set  of  data  files  for  the  periodic 
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orbit,  the  modal  E  matrix,  and  perturbation  series  in  the  form  of  Fourier  series.  The  periodic 
orbit  has  two  large  angles,  Q\  and  02. and  four  small  parameters,  >’2,  V4.  >’5,  and  y$.  The 
change  in  the  two  angles  is  quite  large,  but  the  small  parameters  are  expected  to  be  small  all 
the  time.  The  theory  includes  two  parts.  The  first  part  is  a  data  package  of  Fourier  series, 
which  need  occasional  update.  The  second  part  is  a  set  of  seven  parameters,  <2i,  Q2,  y2,  J4. 
>’5,  y6,  and  B* ,  at  epoch  ta  which  require  more  frequent  update.  Multiply- structured  Fourier 
series  define  the  trajectory  under  modeled  perturbations,  and  it  is  specialized  to  a  particular 
satellite  because  it  is  a  numerical  method.  However,  it  takes  less  than  a  minute  on  a  modern 
personal  home  computer  [69].  After  the  theory  files  are  built  and  the  SGP4  and  TLE 
propagation  is  obtained,  the  low  KAM  torus  is  fitted  to  the  SGP4  and  TLE  data  because 
the  point  on  the  torus  that  matches  the  epoch  time  of  SGP4  and  TLE  data  must  be  found. 
Moreover,  the  fitting  process  corrects  the  frequencies  of  the  torus,  (O \  and  ah,  and  it  figures 
out  the  B*  value  for  the  propagation  period.  It  should  be  noted  that  the  actual  accuracy 
that  can  be  achieved  by  this  theory  is  supposed  to  be  far  better  with  the  raw  observational 
data  than  SGP4  and  TLE  data  due  to  the  inaccuracies  in  the  SGP4  and  TLE.  The  whole 
process  from  the  SGP4  and  TLE  propagation  to  the  least  squares  fitting  needs  no  human 
intervention  because  all  scripts  are  linked  by  bash  files.  The  author  modified  every  script 
to  be  both  compatible  for  Windows™,  and  GNU/Linux  platforms.  The  computational  part 
of  this  effort  has  been  conducted  on  a  64-bit  GNU/Linux  machine. 
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IV.  Results  and  Analysis 


In  this  chapter,  the  results  will  be  presented  and  discussed.  1500  TLEs  were  obtained 
from  www .  space-track .  com  servers,  and  each  were  propagated  for  an  orbit  for  2  months. 
Then,  the  low  KAM  torus  theory  files  were  created,  which  are  associated  with  the  periodic 
orbit,  sectoral  and  tesseral,  eccentricity,  and  air  drag  perturbations.  Finally,  the  low  eccen¬ 
tricity  KAM  torus  was  fitted  to  the  SGP4  and  TLE  propagation  data,  and  the  residuals  were 
calculated.  Therefore,  the  least  squares  fitting  routine  not  only  locates  the  epoch  point  on 
the  KAM  torus  and  corrects  the  frequencies,  but  also  yields  residuals.  The  first  section 
presents  the  success  rate  for  convergence  between  the  low  eccentricity  KAM  torus  and  the 
SGP4  and  TLE  data,  and  provides  details  related  to  the  root  causes  of  failures.  Each  failure 
was  analyzed  individually.  The  second  section  presents  the  relationship  between  orbital 
characteristics  and  residuals  for  all  1500  test  cases.  Again,  some  plots  for  failed  cases  are 
provided  to  figure  out  the  root  causes  of  the  failures.  The  third  section  discusses  and  ana¬ 
lyzes  position  and  velocity  residual  plots  for  LEO  and  GEO  objects  in  the  RAN  reference 
frame,  see  Section  2.4.  The  fourth  section  describes  the  optimal  mean  orbital  limitations 
that  never  fail  the  convergence  of  the  low  eccentricity  KAM  torus.  The  fifth  section  com¬ 
pares  the  low  eccentricity  KAM  torus  and  the  SGP4  predictions  for  the  best  fit  and  a  mean 
test  case.  The  last  section  summarizes  the  chapter. 

4.1  Overview  of  Results 

All  publicly  available  TLEs  which  have  eccentricity  on  the  order  of  10  3  and  smaller 
were  retrieved  from  the  www .  spacetrack .  com  servers.  1500  TLEs  were  pseudo-randomly 
selected  to  fit  the  conditions  mentioned  in  Section  3.1.  First,  each  TLE  was  propagated  for 
an  orbit  for  2  months  period  of  time.  Next,  the  periodic  orbit,  and  its  linearization,  which 
is  Floquet  solution,  for  each  test  case  was  calculated,  and  perturbations  were  added  to  the 
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Table  11.  Rms  Values  for  Least  Squares  Fitting  of  1500  Test  Cases 


Rms  range  (Km) 

Percentage 

0-0.5 

0.6% 

0.5-1 

8.53333% 

1-3 

8.13333% 

3-5 

3.73333% 

5-10 

5.53333% 

10-20 

8.93333% 

20-40 

8.06667% 

40-100 

4.86667% 

100-200 

4.2% 

>200 

4.26667% 

solution.  Then,  the  SGP4  and  TLE  data  were  fit  using  least  squares,  and  the  residuals  for 
each  test  case  were  calculated.  The  success  rate  for  convergence  is  56.8667%.  Table  11 
represents  the  percentage  of  different  range  of  rms  values  for  all  1500  test  cases.  It  should 
be  noted  that  these  percentages  will  increase  if  the  propagation  time  period  is  smaller,  be¬ 
cause  air  drag  is  stochastic,  and  2  months  is  long  enough  to  terminate  the  torus.  Therefore, 
another  analysis  was  conducted.  One  hundred  failed  test  cases  were  pseudo-randomly  se¬ 
lected  and  propagated  for  20  days  instead  of  2  months.  The  success  rate  of  least  squares 
fitting  over  20  days  is  66%.  Table  12  shows  the  rms  values  for  the  100  previously  failed  test 
cases.  It  is  clear  that  the  air  drag  is  the  root  cause  of  the  failure  for  66  test  cases.  The  other 
34  test  cases  will  be  analyzed  individually  in  order  to  find  the  root  causes  of  the  failure. 
Table  13  represents  the  34  failed  test  cases  for  the  second  time. 

Iridium  68,  NORAD  id  25291,  is  an  operational  satellite.  Therefore,  maneuvers  for 
orbit  maintenance,  or  attitude  control  maneuvers  might  cause  the  failure.  Another  reason 
for  the  failure  might  be  its  relatively  low  height  of  semi-major  axis  and  high  B*  values, 
which  are  7. 1558000c +  03  km  and  2.0938000c  —  04,  respectively.  Iridium  68  was  propa¬ 
gated  from  02-01-2013  to  02-09-2013  in  order  to  check  whether  the  air  drag  effect  is  the 
reason  for  the  failure.  However,  the  fitting  process  yields  residuals  on  the  order  of  107  for 


84 


Table  12.  Rms  Values  for  Least  Squares  Fitting  of  the  100  Previously  Failed  Test  Cases 


Rms  range  (Km) 

Percentage 

0-0.5 

0% 

0.5-1 

3% 

1-3 

15% 

3-5 

5% 

5-10 

9% 

10-20 

3% 

20-40 

5% 

40-100 

4% 

100-200 

8% 

>200 

14% 

the  second  time,  which  is  a  typical  error  characteristic  for  polar  inclination.  Therefore,  it 
is  concluded  that  its  inclination  of  8.6397400c  +  01°  causes  the  failure.  In  addition,  all  the 
Iridium  satellites  are  in  a  Walker  constellation,  and  they  have  nearly  the  same  inclinations. 
However,  the  author  tested  all  other  Iridium  debris  and  satellites  to  confirm  whether  the 
high  inclination  is  the  reason  for  the  failure. 

Iridium  33,  NORAD  id  35078,  is  a  debris.  Its  inclination,  semi-major  axis,  and  B* 
values  are  very  close  to  what  Iridium  68,  NORAD  id  25291,  has.  Iridium  33,  NORAD  id 
35078,  was  propagated  from  03-01-2013  to  03-08-2013  in  order  to  check  whether  the  air 
drag  effect  is  the  reason  for  the  failure.  However,  the  residuals  are  identical  to  the  residuals 
that  Iridium  68  yielded  during  the  fitting  process.  Therefore,  it  has  been  concluded  that 
its  inclination  of  8.6399100c +  01°  causes  the  failure.  Moreover,  there  are  4  other  Iridium 
33  debris  which  have  approximately  the  same  inclinations.  Thus,  all  Iridium  33  debris  in 
Table  13  fail  due  to  their  high  inclinations. 

Meteor  1-21,  NORAD  id  7714,  is  an  operational  weather  satellite.  Its  period  is  6.1373891c  + 
03  seconds,  which  is  in  the  close  vicinity  of  the  14: 1  resonance.  It  was  propagated  from  09- 
02-2013  to  09-08-2013.  It  yielded  residuals  on  the  order  of  hundreds  of  meters  during  the 
fitting  process  although  the  convergence  was  accomplished.  Because  the  fit  was  converged, 
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Table  13.  34  Test  Cases  Failed  for  the  Second  Time 


NORAD  id 

Satellite  Name 

Type 

Launch  Date 

25291 

IRIDIUM  68 

PAYLOAD 

1998-04-07 

35078 

IRIDIUM  33 

DEBRIS 

1997-09-14 

7714 

METEOR  1-21 

PAYLOAD 

1975-04-01 

29349 

KOREASAT  5 

PAYLOAD 

2006-08-22 

25531 

IRIDIUM  83 

PAYLOAD 

1998-11-06 

26599 

BEIDOU  1 

PAYLOAD 

2000-10-30 

34356 

IRIDIUM  33 

DEBRIS 

1997-09-14 

16495 

COSMOS  1726 

PAYLOAD 

1986-01-17 

33697 

FENGYUN  1C 

DEBRIS 

1999-05-10 

9701 

DELTA  1 

DEBRIS 

1973-11-06 

21153 

SL-8  R/B 

ROCKET  BODY 

1991-03-12 

24115 

PEGASUS 

DEBRIS 

1994-05-19 

23092 

COSMOS  2279 

PAYLOAD 

1994-04-26 

31039 

FENGYUN  1C 

DEBRIS 

1999-05-10 

15308 

GALAXY  3 

PAYLOAD 

1984-09-21 

28923 

FREGAT  R/B 

ROCKET  BODY 

1991-03-12 

31988 

FENGYUN  1C 

DEBRIS 

1999-05-10 

34928 

IRIDIUM  33 

DEBRIS 

1997-09-14 

30760 

FENGYUN  1C 

DEBRIS 

1999-05-10 

31561 

FENGYUN 1C 

DEBRIS 

1999-05-10 

29830 

FENGYUN  1C 

DEBRIS 

1999-05-10 

31339 

FENGYUN 1C 

DEBRIS 

1999-05-10 

33605 

FENGYUN  1C 

DEBRIS 

1999-05-10 

29921 

FENGYUN 1C 

DEBRIS 

1999-05-10 

20061 

NAVSTAR  14 

PAYLOAD 

1989-06-10 

7560 

SL-8  DEB 

DEBRIS 

1972-08-16 

36482 

IRIDIUM  33 

DEBRIS 

1997-09-14 

13073 

COSMOS  1275 

DEBRIS 

1981-06-04 

24435 

EXPRESS  2 

PAYLOAD 

1996-09-26 

21782 

COSMOS  2168 

PAYLOAD 

1991-11-12 

14966 

SL-8  R/B 

ROCKET  BODY 

1984-05-11 

35981 

COSMOS  2251 

DEBRIS 

1993-06-16 

34488 

IRIDIUM  33 

DEBRIS 

1997-09-14 

23404 

COSMOS  2297 

PAYLOAD 

1994-11-24 
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it  was  propagated  for  some  random  dates  over  5  days  in  order  to  observe  an  improvement 
in  the  residuals.  However,  the  residuals  appeared  to  be  poor  for  all  attempts.  Its  semi-major 
axis  and  B*  values  can’t  cause  such  poor  fits.  Therefore,  it  is  concluded  that  the  resonance 
is  the  reason  for  the  poor  residuals.  The  KAM  torus  might  deform  rapidly  near  the  resonant 
orbits  if  not  fail  completely. 

Koreasat  5,  NORAD  id  29349,  is  an  operational  satellite.  It  was  propagated  from  05- 
01-2013  to  05-10-2013  and  from  07-01-2013  to  07-08-2013  in  order  to  check  whether  the 
third  body  mass  perturbations  lead  to  the  failure.  Both  propagations  yielded  the  same  big 
residuals  during  the  fitting  process.  Its  period  is  8.6165984c  +  04  seconds,  and  it  is  in  the 
close  vicinity  of  the  1:1  resonance.  Therefore,  it  is  concluded  that  the  resonance  is  the 
reason  for  the  failure. 

Iridium  83,  NORAD  id  2553 1 ,  is  an  operational  satellite.  It  was  propagated  from  02-01- 
2013  to  02-08-2013.  The  fitting  process  yielded  characteristic  polar  inclination  residuals, 
which  were  on  the  order  of  107  km.  Therefore,  it  is  concluded  that  high  inclination,  which 
is  8.6396500c  +  01°,  is  the  reason  for  the  failure. 

Beidou  1,  NORAD  id  26599,  is  an  operational  satellite.  It  was  propagated  from  01-01- 
2013  to  01-10-2013  and  from  09-01-2013  to  09-10-2013.  Both  least  squares  fits  yielded 
big  residuals,  and  the  convergence  failed  for  both  attempts.  The  resonance  might  be  the 
reason  for  the  failure  because  its  period  is  8.7352213c  +  04  seconds. 

Cosmos  1726,  NORAD  id  16495,  is  an  operational  satellite.  It  was  propagated  from 
01-01-2014  to  01-14-2014.  It  has  an  altitude  of  527  km,  and  aB*  value  of  3.4805000c  — 04. 
The  orbit  was  successfully  fitted  to  SGP4  and  TLE  with  an  rms  value  of  7.54832  km  over 
13  days.  The  poor  fit  shows  that  the  atmospheric  model  used  in  the  program,  which  is 
based  on  the  U.S.  Standard  Atmosphere,  should  be  improved  to  predict  satellites  with  low 
altitudes  with  better  accuracy.  Figure  23  represents  the  plots  for  position  and  velocity 
residuals.  In  the  plot,  along-track  residuals  are  relatively  bigger  than  the  other  components, 
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Figure  23.  The  Residuals  for  Cosmos  1726  over  13  days 


which  proves  that  the  air  drag  effect  is  the  reason  for  the  failure  of  the  first  attempt. 

Fengyun  1C  debris,  NORAD  id  33697,  31988,  30760,  31561,  31339,  33605,  and 
29921,  are  in  the  close  vicinity  of  the  14:1  resonance.  Each  debris  was  propagated  for 
a  short  period  of  time  in  order  to  figure  out  whether  the  air  drag  effect  is  the  reason  for  the 
failure.  None  of  the  fits  converged.  Therefore,  the  reason  for  the  failure  is  the  14:1  reso¬ 
nance.  Moreover,  Fengyun  1C  debris,  NORAD  id  31039,  and  29830,  prove  this  statement 
because  although  two  of  them  have  lower  height  and  higher  B *  value,  both  fits  converged, 
and  the  residuals  weren’t  bad  for  their  low  altitude.  Fengyun  1C  debris,  NORAD  id  29830, 
yields  an  rms  value  of  9.47329  km,  and  NORAD  id  31039  yields  an  rms  value  of  22.8057 
km.  Figure  24  represents  the  plot  for  position  and  velocity  residuals  of  Fengyun  1C  de¬ 
bris,  NORAD  id  29830,  and  Figure  25  shows  the  plot  for  position  and  velocity  residuals  of 
Fengyun  1C  debris,  NORAD  id  31039.  In  the  plots,  the  quadratic  characteristic  of  air  drag 
effect  can  be  easily  seen. 

Delta  1  debris,  NORAD  id  9701,  has  relatively  high  B*  value,  which  is  1. 2101000c  — 
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Figure  24.  The  Residuals  for  Fengyun  1C,  NORAD  id  29830,  over  6  Days 
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Figure  25.  The  Residuals  for  Fengyun  1C,  NORAD  id  31039,  over  6  Days 
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02,  and  it  is  also  in  the  close  vicinity  of  the  13:1  resonance.  It  was  propagated  from  01- 
01-2014  to  01-09-2014.  The  least  squares  fit  converged  with  an  rms  value  of  94  km  over  8 
days. Therefore,  the  air  drag  effect  is  believed  to  be  the  prime  reason  for  the  failure. 

SL-8  R/B,  NORAD  id  21153,  is  a  rocket  body.  It  has  an  inclination  of  8.2926500c  + 
01°,  a  semi-major  axis  of  7.3573300e  +  03  km,  a  B*  value  of  1.3800000c  — 04,  and  a  period 
of  6.2804555c +  03  seconds.  It  was  propagated  from  01-02-2014  to  01-08-2014.  The  least 
squares  fitting  process  converged  with  an  rms  value  of  47.0068  km.  Thus,  the  air  drag  is 
the  reason  for  the  failure.  Figure  26  represents  the  plot  for  position  and  velocity  residuals. 
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Figure  26.  The  Residuals  for  SL-8  R/B,  NORAD  id  21153,  over  10  Days 


Pegasus,  NORAD  id  24115,  is  a  debris.  Its  relatively  low  altitude,  which  is  512.275 
km,  and  relatively  high  B*  value,  which  is  6.3281000c  —  04,  are  responsible  for  the  failure 
because  the  least  squares  fit  converged  when  it  was  propagated  from  01-02-2014  to  01-09- 
2014  over  8  days.  Figure  27  represents  the  plot  for  position  and  velocity  residuals. 

Cosmos  2279,  NORAD  id  23092,  is  an  operational  satellite.  It  was  propagated  from 
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Figure  27.  The  Residuals  for  Pegasus,  NORAD  id  24115,  over  9  Days 


01-01-2014  to  01-19-2014.  The  least  squares  fit  converged  with  an  rms  value  of  13.8748 
km  over  18  days.  Therefore,  it  is  clear  that  the  least  squares  fit  has  failed  to  converge  due 
to  the  air  drag  effect.  Figure  28  represents  the  plot  for  position  and  velocity  residuals. 

Galaxy  3,  NORAD  id  15308,  is  an  operational  satellite  with  a  period  of  24.03.  The 
propagation  from  01-01-2014  to  01-10-2014  yielded  an  rms  value  of  1.94622  km.  There¬ 
fore,  the  inaccuracies  in  the  TLEs  caused  the  failure. 

Fregat  R/B,  NORAD  id  28923,  is  a  rocket  body.  It  has  an  inclination  of  5.6419 lOOe  + 
01°,  a  semi-major  axis  of  2.9813700^  +  04  km,  and  a  period  of  5.1231209e  +  04  seconds. 
It  was  propagated  for  different  time  periods,  and  all  attempts  yielded  characteristic  rms 
values  for  resonant  orbits,  which  were  approximately  on  the  order  of  10,  and  the  residuals 
were  random-like.  Therefore,  it  is  concluded  that  the  resonance  is  the  reason  for  the  failure. 

Navstar  14,  NORAD  id  20061,  is  an  operational  satellite.  It  has  an  inclination  of 
5.5107800^  +  01°,  a  semi-major  axis  of  2.7647000^  +  04  km,  and  a  period  of  4.5749 126e  + 
04  seconds.  It  was  propagated  from  01-01-2014  to  01-10-2014.  The  least  squares  fitting 
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Figure  28.  The  Residuals  for  Cosmos  2279,  NORAD  id  23092,  over  18  Days 


process  yielded  random  big  rms  values  for  iterations,  which  were  on  the  order  of  10  6 . 
Therefore,  it  is  concluded  that  the  resonance  is  the  reason  for  the  failure. 

SL-8  debris,  NORAD  id  7560,  has  a  period  of  6. 1487649c  +  03  seconds.  Thus,  it  is  in 
the  close  vicinity  of  the  14:1  resonance.  It  was  propagated  from  01-01-2014  to  01-10-2014. 
The  least  squares  fit  failed  to  converge,  and  yielded  random  big  rms  values  for  iterations, 
which  is  typical  residuals  for  orbits  in  the  close  vicinity  of  the  resonance. 

Cosmos  1275,  NORAD  id  13073,  is  a  debris.  It  was  propagated  from  01-02-2014  to  01- 
09-2014.  The  convergence  was  accomplished  with  an  rms  value  of  3.39862  km.  Therefore, 
the  reason  for  the  failure  is  the  air  drag  effect.  Figure  29  represents  the  plots  for  position 
and  velocity  residuals  over  7  days. 

Express  2,  NORAD  id  24435,  is  an  operational  satellite.  It  has  an  inclination  of 
1.2856800e  +  01°,  a  semi-major  axis  of  4. 2 182300c  +  04  km,  and  a  period  of  8. 6219633c  + 
04  seconds.  It  is  in  the  close  vicinity  of  the  1:1  resonance.  Express  2  was  propagated  from 
01-03-2014  to  01-08-2014  over  5  days.  It  yielded  random  big  rms  values  during  the  least 
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Figure  29.  The  Residuals  for  Cosmos  1275,  NORAD  id  13073,  over  7  Days 

squares  fitting  process,  and  it  failed  to  converge.  Therefore,  the  resonance  is  the  reason  for 
the  failure  not  the  third  body  mass  perturbations. 

Cosmos  2168,  NORAD  id  21782,  is  an  operational  satellite.  It  has  an  inclination  of 
8.2598100e  +  01°,  a  semi-major  axis  of  7.7788100^  +  03  km,  a  B*  value  of  2.4820000e  — 
04,  and  a  period  of  6.8277959^  +  03  seconds.  It  was  propagated  from  01-02-2014  to  01- 
09-2014  over  7  days.  The  least  squares  fit  converged  with  an  rms  value  of  2.11935  km. 
Therefore,  it  is  concluded  that  the  reason  for  the  failure  is  the  air  drag  effect.  Figure  30 
represents  the  plots  for  position  and  velocity  residuals  over  7  days. 

SL-8  R/B,  NORAD  id  14966,  has  an  inclination  of  8.2972500e  +  01°,  a  semi-major  axis 
of  7.3620200e  +  03  km,  a  B*  value  of  1.9141000e  —  04,  and  a  period  of  6.2864617e  +  03 
seconds.  It  was  propagated  from  01-01-2014  to  01-09-2014  over  9  days.  The  least  squares 
fit  converged  with  an  rms  value  of  13.9352  km.  Therefore,  the  air  drag  effect  is  the  reason 
for  the  failure. 

Cosmos  2251,  NORAD  id  35981,  has  an  inclination  of  7.4059500e  +  01°,  a  semi-major 
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Figure  30.  The  Residuals  for  Cosmos  2168,  NORAD  id  21782,  over  9  Days 


axis  of  6.9195900e  +  03  km,  a  B*  value  of  2.9060000e  —  03,  and  a  period  of  5.7283735e  + 
03  seconds.  It  was  propagated  from  01-01-2013  to  01-14-2013  over  13  days.  The  least 
squares  fit  converged  with  an  rms  value  of  86.645  km,  which  isn’t  surprising  because  of  its 
low  altitude  and  high  B*  value.  Therefore,  the  air  drag  effect  is  the  reason  for  the  failure. 
Figure  31  represents  the  plots  for  position  and  velocity  residuals  over  15  days. 

Cosmos  2297,  NORAD  id  23404,  is  in  the  close  vicinity  of  the  14:1  resonance.  It  has 
an  inclination  of  7. 1013200c  +  01°,  a  semi-major  axis  of  7.2261100c  +  03  km,  a  B*  value 
of  2.4820000c  —  04,  and  a  period  of  6. 113 1864c +  03  seconds.  It  was  propagated  from 
01-01-2014  to  01-15-2014,  01-09-2014,  and  01-07-2014  for  three  different  time  periods. 
All  attempts  failed  because  of  the  resonance. 

In  conclusion,  1500  test  cases,  which  fit  the  conditions  mentioned  in  Section  3.1,  were 
propagated  using  the  low  eccentricity  KAM  torus  method.  Then,  each  KAM  torus  propa¬ 
gation  data  was  fitted  to  SGP4  and  TLE  prediction  data  associated  with  each  test  case  over 
2  months.  Initial  success  rate  for  convergence  is  56.8667%.  Then,  100  new  test  cases  were 
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Figure  31.  The  Residuals  for  Cosmos  2251,  NORAD  id  35981,  over  15  Days 


selected  from  the  failed  test  cases,  and  the  same  process  was  applied  to  100  new  test  cases. 
The  success  rate  of  least  squares  fitting  over  20  days  instead  of  2  months  for  the  new  test 
cases  is  66%.  Next,  34  failed  test  cases,  which  failed  to  converge  for  the  second  time,  were 
individually  analyzed.  Seven  test  cases  failed  for  the  third  time  because  their  inclinations 
are  in  the  close  vicinity  of  polar  inclination.  An  inclination  of  86°  caused  the  construction 
of  the  KAM  torus  to  fail.  The  unstable  forms  of  the  Legendre  polynomial  recursions  for 
the  calculation  of  geopotential  caused  the  high  eccentricity  problem.  The  solution  is  left 
for  future  studies.  One  test  case  failed  to  converge  to  the  SGP4  and  TLE  data  because  of 
the  inaccuracies  in  the  TLEs.  Ten  test  cases  failed  due  to  the  air  drag  effect.  Test  cases  with 
relatively  high  B*  value  and  low  altitude  can  only  be  propagated  for  short  period  of  time 
because  of  effective  air  drag.  Sixteen  test  cases  failed  to  converge  due  to  the  resonance. 
There  is  no  KAM  torus  near  the  resonance.  Therefore,  a  clever  approach  should  be  taken 
to  deal  with  the  resonance,  which  is  important  for  the  generalization  of  the  theory  to  all 
orbits.  In  addition,  critical  and  polar  inclination  problems  should  also  be  addressed  for  the 
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generalization  of  the  theory. 


4.2  Orbital  Characteristics  and  Residuals 

The  relationship  between  some  orbital  characteristics,  such  as  semi-major  axis,  eccen¬ 
tricity,  inclination,  B*  value,  and  period,  and  rms  values  yield  valuable  information  about 
the  theory.  Because  many  factors  degrade  the  accuracy  and  the  dimensions  that  can  be 
used  to  present  the  data  are  restricted  to  3,  straightforward  presentation  of  data  isn’t  useful. 
Therefore,  success  ratio,  which  is  successful  attempts  over  the  sum  of  the  successful  and 
failed  attempts,  for  fixed  time  intervals  versus  orbital  characteristics  are  used  for  the  plots. 
The  time  intervals  vary  based  on  the  orbital  characteristics.  However,  there  are  30  time 
intervals  for  all  orbital  characteristics.  Each  success  ratio  is  an  average  value  for  a  time  in¬ 
terval,  which  is  attributed  to  the  starting  point  of  the  interval.  This  is  important  information 
for  evaluating  the  plots  including  resonance  data  because  the  effects  of  the  resonance  don’t 
seem  to  be  aligned.  Therefore,  additional  explanations  are  provided  for  the  plots.  The  data 
obtained  from  1500  test  cases  were  used  for  all  plots  in  Section  4.1. 

The  success  ratio  versus  semi-major  axis  was  plotted  for  LEO  and  GEO  satellites  be¬ 
cause  there  were  few  test  cases  at  the  MEO  altitudes.  For  LEO  satellites,  the  prime  pertur¬ 
bation  that  degrades  the  accuracy  for  the  theory  is  the  air  drag.  The  other  parameter  that 
degrades  the  accuracy  for  LEO  satellites  is  the  resonance.  For  GEO  satellites,  the  prime 
perturbation  is  the  third  body  mass.  The  other  parameter  that  degrades  the  accuracy  for 
GEO  satellites  is  the  1:1  resonance. 

Figure  32  represents  the  success  ratio  versus  semi-major  axis  for  the  test  cases  which 
have  semi-major  axes  of  8000  km  and  smaller.  The  resonance  points  from  left  to  right  are 
the  29:2,  14:1,  27:2,  13:1,  37:3,  and  12:1  resonances.  Although  29:2  and  37:3  resonances 
do  seem  to  be  aligned  with  the  decrease  in  the  success  ratio,  the  success  ratio  is  an  average 
over  a  time  interval  attributed  to  the  starting  point  of  the  interval.  For  the  29:2  resonance 
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case,  the  success  ratio  average  is  from  7050  km  to  7100  km,  and  the  altitude,  which  is 
calculated  from  the  center  of  Earth,  for  the  29:2  resonance  is  7090.846  km  for  orbits  with 
zero  eccentricities.  For  the  37:3  resonance  case,  the  success  ratio  average  is  from  7850 
km  to  7900  km,  and  the  altitude,  which  is  calculated  from  the  center  of  Earth,  for  the  37:3 
resonance  is  7898.714  km  for  orbits  with  zero  eccentricities.  In  addition,  the  13:1  resonance 
doesn’t  degrade  the  accuracy  to  zero  for  the  interval  between  7600  km  and  7700  km.  There 
are  missing  data  for  the  interval  between  7600  km  and  7650  km.  However,  the  resonance 
certainly  degrades  the  accuracy.  The  success  ratio  increases  with  increasing  altitude,  which 
isn’t  surprising,  because  air  drag  perturbation  is  more  effective  at  low  altitudes. 
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Figure  32.  The  Success  Ratio  versus  Semi-Major  Axis  for  LEO  Objects  over  2  Months 

The  plot  for  the  success  ratio  of  GEO  objects  is  surprising  at  the  first  glimpse  because 
there  are  unexpected  failures.  Further  examination  yielded  that  three  test  cases  created 
these  failures.  The  author  propagated  these  three  test  cases,  which  are  Ius  R/B,  NORAD 
id  22316,  Tvsat  2,  NORAD  id  20168,  and  familiar  Galaxy  3,  NORAD  id  15308,  from 
01-01-2014  to  03-01-2014.  They  converged  with  rms  values  of  13.3379  km,  10.1593  km, 
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and  14.4785  km,  respectively,  which  proves  that  many  failures  can  be  related  to  the  inac¬ 
curacies  in  TLEs,  but  this  effort  ignores  this  fact.  Figure  33  represents  the  success  ratio  for 
test  cases  which  have  a  semi-major  axis  between  42100  km  and  42500  km  before  and  after 
corrections  are  applied.  The  accuracy  decreases  rapidly  after  42316  km.  The  convergence 
totally  vanishes,  or  the  accuracy  becomes  extremely  low  if  the  convergence  is  accomplished 
for  42416  km  and  bigger  semi-major  axis  values  due  to  the  third  body  mass  perturbations. 
The  1:1  resonance  certainly  degrades  the  accuracy.  In  addition,  active  payloads  in  geosyn¬ 
chronous  orbit  are  maneuvered  frequently  to  suppress  the  action  of  the  resonance. 
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Figure  33.  The  Success  Ratio  versus  Semi-Major  Axis  for  GEO  Objects 


Figure  34  represents  the  success  ratio  versus  eccentricity  for  the  test  cases  which  have 
eccentricities  between  0  and  0.01.  In  the  plot,  from  0  to  0.004  the  accuracy  is  pretty  lev¬ 
eled,  except  some  ups  and  downs  due  to  other  orbital  characteristics.  However,  after  the 
threshold  of  0.004  eccentricity  is  passed  the  success  ratio  decreases  to  approximately  30%. 
This  proves  how  eccentricity  is  an  important  factor  in  terms  of  accuracy  for  the  theory. 
Moreover,  another  analysis  in  Section  4.4  will  examine  the  upper  limit  for  the  eccentricity 
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that  can  be  modeled  by  the  theory. 


1 

0.8 

o 

S  0-6 

DC 

<n 

</> 

0) 

o 

3  0.4 

( I, ) 

0.2 

0 

0  0.001  0.002  0.003  0.004  0.005  0.006  0.007  0.008  0.009  0.01 

Eccentricity 

Figure  34.  The  Success  Ratio  versus  Eccentricity  over  2  Months 

Figure  35  represents  the  success  ratio  versus  B*  for  the  test  cases  which  have  B*  values 
between  0  and  0.01.  In  the  plot,  when  the  B*  value  of  0.003  is  passed  the  success  ratio 
for  convergence  decreases  to  approximately  zero.  The  sharp  decrease  with  increasing  B* 
value  shows  the  importance  of  the  value  for  the  theory.  As  mentioned  previously,  the  B* 
value  defines  the  susceptibility  of  near-Earth  objects  to  air  drag.  The  more  effective  the 
air  drag  is  on  a  near-Earth  object,  the  less  accurate  low  eccentricity  KAM  torus  theory 
becomes.  The  drastic  effect  of  the  B*  value  on  the  theory  necessitates  either  an  improved 
atmospheric  model  or  the  shorter  propagation  time  interval.  However,  the  first  option  may 
not  introduce  enough  improvement  in  the  accuracy  due  to  the  stochastic  nature  of  air  drag. 

One  thousand  five  hundred  test  cases  were  pseudo-randomly  selected  under  the  conditions 
mentioned  in  details  in  Section  3.1,  which  include  the  close  vicinity  of  polar  and  critical 
inclination,  and  eccentricity  on  the  order  of  10  3  and  smaller.  Figure  36  represents  the 
success  ratio  versus  inclination  in  order  to  check  whether  there  are  other  inclination  values 
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Figure  35.  The  Success  Ratio  versus  Bstar  over  2  Months 
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Figure  36.  The  Success  Ratio  versus  Inclination  over  2  Months 
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that  lead  to  the  resonance.  In  the  plot,  the  inclinations  of  20°  and  80°  seem  to  decrease  the 
accuracy  drastically,  but  this  is  an  unexpected  situation.  Further  examination  yielded  that 
three  test  cases  for  inclinations  between  20°  and  25°,  two  test  cases  for  inclinations  between 
25°  and  30°,  and  three  test  cases  for  inclinations  between  75°  and  80°  have  created  these 
failures.  For  inclinations  20°-25°,  Delta  2  rocket  body,  NORAD  id  21931,  Pegasus  rocket 
body,  NORAD  id  22491,  and  Scd  2,  NORAD  id  25504,  have  caused  the  failures.  They  all 
have  low  altitudes.  Delta  2  R/B  was  propagated  from  01-01-2013  to  01-10-2013,  and  the 
other  two  were  propagated  from  01-01-2014  to  01-10-2014.  They  converged  to  the  SGP4 
and  TLE  data  with  rms  values  of  15.2938  km,  0.3447  km,  and  0.2793  km,  respectively. 
For  inclinations  25°-30°,  Tansei  1,  NORAD  id  4952,  Delta  1  debris,  NORAD  id  10309, 
and  Galex,  NORAD  id  27783,  caused  the  failures.  Delta  1  and  Galex  have  low  altitudes. 
Two  of  them  were  propagated  from  01-01-2014  to  01-10-2014  with  rms  values  of  269.943 
km,  0.2481  km,  respectively.  Although  Delta  1  successfully  converged,  its  semi-major 
axis,  which  is  6849.32,  caused  a  big  residual.  Tansei  1  was  also  propagated  from  from 
01-01-2014  to  01-10-2014  in  order  to  check  whether  attitude  control,  or  station  keeping 
maneuvers  caused  the  failure.  It  was  converged  with  an  rms  value  of  0.5388  km.  For 
inclinations  75°-80°,  Scout  X-4,  NORAD  id  871,  and  SL-21  rocket  body,  NORAD  id 
29159,  caused  the  failures.  Two  of  them  were  propagated  from  01-01-2014  to  01-10-2014 
with  rms  values  of  0.9033  km,  and  24.7307  km,  respectively.  Therefore,  the  polar  and 
the  critical  inclinations  are  the  only  ones  that  fail  the  converge  and  degrade  accuracy.  The 
success  ratio  is  quite  leveled  if  the  corrections  are  applied.  Thus,  there  is  no  linear  relation 
between  the  success  ratio  and  the  inclination.  Moreover,  another  analysis  in  Section  4.4 
will  examine  the  limits  for  the  polar  and  critical  inclinations  that  can  be  modeled  by  the 
theory. 

Figure  37  represents  the  success  ratio  versus  the  period  for  the  test  cases  which  have 
periods  between  6000  seconds  and  7160  seconds.  The  resonance  points  from  left  to  right 
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are  the  14:1,  13:1,  and  12:1.  The  other  order  resonances,  such  as  27:2  and  37:3,  can’t  be 
detected  in  the  plot.  The  two  drastic  decreases  between  6400  seconds  and  6600  seconds  are 
due  to  the  small  number  of  test  cases  which  all  except  one  have  high  B*  values.  One  of  them 
failed  to  converge  due  to  the  inaccurate  TLEs.  For  periods  between  6440  seconds  and  6480 
seconds.  Delta  1  rocket  body,  NORAD  id  21602,  and  Nimbus  2  debris,  NORAD  id  31579, 
reduced  the  accuracy.  For  periods  between  6520  seconds  and  6560  seconds,  Cosmos  770, 
NORAD  id  8325,  which  has  a  B*  value  of  2.845  x  10-3,  Meteor  3-2,  NORAD  is  19336, 
which  has  a  B*  value  of  1.1514  x  10^3,  Cosmos  1275  debris,  NORAD  id  14440,  which  has 
a  B*  value  of  1.4665  x  10-2,  and  SF-8  rocket  body,  NORAD  id  11170,  which  failed  due  to 
the  inaccurate  TFEs  because  propagation  from  01-01-2002  to  03-01-2002  yielded  16.7697 
km  residual.  Because  there  are  few  test  cases  between  6400  seconds  and  6600  seconds,  the 
resonance  effect  can’t  be  detected  during  this  time  period.  Therefore,  if  the  propagation 
time  is  decreased,  the  success  ratio  between  6400  seconds  and  6600  seconds  has  a  very 
similar  trend  as  the  success  ratio  between  6200  seconds  and  6400  seconds,  that  is  a  linear 
increase  with  smaller  variations. 

In  conclusion,  the  resonance  has  emerged  as  a  problem  again  in  the  broader  picture. 
This  section  presents  the  analysis  for  1500  test  cases.  For  FEO  objects,  air  drag  degrades 
the  success  ratio  of  the  least  squares  fit  over  2  months  to  zero  for  altitudes  below  422 
km.  The  accuracy  also  increases  with  increasing  altitude.  For  GEO  satellites,  third  body 
mass  perturbation  degrades  the  success  ratio  of  the  fit  to  zero  for  altitudes  above  36038 
km.  The  only  resonance  that  effects  the  accuracy  is  the  1:1  resonance  for  altitudes  between 
35622  km  and  36122  km.  The  eccentricity  has  appeared  as  an  important  parameter  that 
degrades  the  accuracy.  A  sharp  decrease  in  the  accuracy  occurs  for  eccentricities  bigger 
than  0.004.  The  B*  value  has  drastically  decreased  the  accuracy  of  the  theory.  Moreover, 
all  16  test  cases  between  the  interval  0.0030  and  0.0033  have  failed  to  converge.  Therefore, 
the  upper  limit  for  B*  for  the  theory  can  be  defined  as  0.003.  However,  the  success  rate  for 
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Figure  37.  The  Success  Ratio  versus  Period  over  2  Months 


B*  values  between  0  and  0.001  can  be  easily  enhanced  by  an  improved  atmospheric  model. 
For  the  current  theory,  217  out  of  725  test  cases  have  failed  for  the  B*  values  between  0 
and  0.00033.  The  polar  inclination  and  the  critical  inclination  were  proved  to  be  the  only 
problem  that  terminates  the  convergence  of  the  least  squares  fit.  Thus,  there  are  no  other 
inclination  values  that  lead  to  failure  of  the  convergence  of  the  least  squares  fit.  Moreover, 
there  is  no  linear  relation  between  the  success  ratio  and  the  inclination.  The  success  ratio 
based  on  the  period  is  unsurprisingly  very  similar  to  the  success  ratio  based  on  the  altitude. 
There  is  a  linear  relation  between  the  period  and  the  success  ratio  of  the  convergence  of  the 
least  squares  fit. 


4.3  Some  Samples  from  the  Results 

The  general  performance  of  the  theory  has  been  presented  so  far.  This  section  intro¬ 
duces  some  characteristic  examples  of  the  least  squares  fits  for  the  LEO  and  GEO  objects. 
Every  test  case  has  a  plot  that  represents  the  accuracy  of  the  low  eccentricity  KAM  torus 


103 


theory  assuming  that  SGP4  and  TLE  data  are  raw  observational  data.  The  reference  frame 
for  the  plots  is  RAN  frame  of  reference,  see  Section  2.4.  However,  this  isn’t  a  fair  evalu¬ 
ation  of  the  accuracy  of  the  theory  because  the  basis  frequencies  and  B*  are  continuously 
updated  using  SGP4  and  TLE  data.  Section  4.5  introduces  a  fair  comparison  between 
SGP4  and  the  low  eccentricity  KAM  torus  method  in  terms  of  accuracy.  Therefore,  the 
plots  introduced  in  this  section  present  whether  the  numerically  built  torus  can  be  fitted  to 
observational  data.  If  the  torus  is  fitted  to  the  observational  data  with  high  accuracy,  the 
result  is  a  compact  numerical  set  of  algorithms  that  are  as  accurate  as  numerical  methods 
and  as  fast  as  general  perturbation  techniques,  in  terms  of  computational  time. 
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Figure  38.  The  Residuals  for  Thorad  Delta  1  Debris,  NORAD  id  8168,  over  2  Months 


Figure  38  and  Figure  39  represent  the  best  2  least  squares  fits  for  LEO  objects.  They 
are  Thorad  Delta  1  debris,  NORAD  is  8168,  and  Thorad  Delta  1  debris,  NORAD  is  8140. 
The  least  squares  fitting  process  yielded  3. 6864262c  —  01  km  and  4.2330236c  —  01  km, 
respectively.  Although  both  of  them  have  similar  orbital  characteristics,  the  difference 
in  the  accuracy  is  due  to  the  B *  values,  which  are  8.2174000c  —  05  and  9.3318000c  — 
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Figure  39.  The  Residuals  for  Thorad  Delta  1  Debris,  NORAD  id  8140,  over  2  Months 

05,  respectively.  It  should  be  also  noted  that  Thorad  Delta  1  debris,  NORAD  is  8140,  is 
approximately  at  an  altitude  25  km  higher  than  Thorad  Delta  1  debris,  NORAD  is  8168. 
Both  plots  present  the  bounded  oscillatory  behavior  of  the  residuals  over  2  months.  They 
both  have  an  altitude  slightly  above  1450  km.  Thus,  air  drag  isn’t  that  effective,  and  the 
quadratic  structure  of  air  drag  usually  appeared  on  along-track  and  radial  components  can’t 
be  detected  in  the  plots. 

Ops  6630  debris,  NORAD  id  9323,  Cosmos  2251  debris,  NORAD  id  35469,  and  Tho¬ 
rad  Agena  D  debris,  NORAD  id  4214,  have  least  squares  fits  with  residuals  on  the  order 
of  kilometers  due  to  different  reasons.  Figure  40  represents  the  position  and  the  velocity 
residuals  of  Ops  6630  debris  for  the  least  squares  fit.  Ops  6630  debris,  NORAD  id  9323, 
has  an  inclination  of  9.7031200e  +  01,  a  semi-major  axis  of  7.8056200e  +  03,  a  B*  value 
of  2. 5971000c  — 04,  and  a  period  of  6.863 1247c +  03  seconds.  Its  relatively  high  B*  value 
leads  to  a  residual  of  5.02942  km  over  2  months.  However,  the  familiar  quadratic  structure 
of  air  drag  that  appears  on  the  along-track  and  the  radial  components  can’t  be  detected  in 
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Figure  40.  The  Residuals  for  Ops  6630  Debris,  NORAD  id  9323,  over  2  Months 


the  plot  because  its  B*  value  doesn’t  exceed  the  limit  that  can  be  handled  by  the  theory. 
However,  Figure  41  represents  the  quadratic  air  drag  effect  on  the  least  squares  fit  of  Tho- 
rad  Agena  D  debris.  The  along-track  and  the  radial  components  seem  to  diverge  from  the 
desired  oscillatory  behavior  of  the  residuals.  Another  factor  that  degrades  the  accuracy  is 
the  inaccuracies  in  the  TLEs.  Figure  42  represents  the  residuals  of  Cosmos  2251  debris  for 
the  least  squares  fit.  In  the  plot,  one  of  the  TLEs  is  certainly  inaccurate. 

GEO  satellites  don’t  yield  neat  oscillatory  behavior  for  the  residuals  because  the  third 
body  mass  perturbations  aren’t  added  to  the  theory  yet.  However,  third  body  mass  pertur¬ 
bations  can  be  calculated  like  air  drag  perturbations,  see  Section  3.5.  There  will  be  at  least 
one  additional  angle  in  the  Fourier  series  to  specify  the  motion  of  the  Sun  or  the  Moon. 
Three  example  test  cases  were  chosen  to  present  the  characteristic  fits  for  GEO  satellites. 
Two  of  them  have  the  best  fits  for  GEO.  The  other  one  is  one  of  the  poorest  fits.  The  poor 
fit  proves  the  behavior  of  the  residuals  when  the  altitude  of  the  satellite  approaches  to  the 
threshold  of  36000  km.  Figure  43  represents  the  position  and  the  velocity  residuals  of  Ops 
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Figure  41.  The  Residuals  for  Thorad  Agena  D  Debris,  NORAD  id  4214,  over  2  Months 
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Figure  42.  The  Residuals  for  Cosmos  2251  Debris,  NORAD  id  35469,  over  2  Months 
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Figure  43.  The  Residuals  for  Ops  9443  payload,  NORAD  id  11621,  over  2  Months 

9443  payload,  NORAD  id  11621,  for  the  least  squares  fit.  The  least  squares  fit  for  Ops 
9443  payload  was  converged  with  an  rms  value  of  8.4179  km.  Although  this  fit  is  one  of 
the  best  fits  for  GEO  satellites,  its  accuracy  is  far  worse  than  the  best  fit  for  LEO  objects 
because  objects  at  high  altitudes  change  their  orbits  towards  the  ecliptic  plane  of  the  solar 
system.  Therefore,  the  third  body  mass  perturbations  should  be  added  to  the  theory  for  bet¬ 
ter  accuracy  for  GEO  satellites.  Figure  44  represents  the  position  and  the  velocity  residuals 
of  Ops  9438  payload,  NORAD  id  10001,  for  the  least  squares  fit.  There  are  irregular  pat¬ 
terns  in  the  plot.  However,  third  body  mass  perturbations  are  not  the  only  source  of  error. 
Both  station  keeping  maneuvers  and  inaccurate  TLEs  can  cause  these  irregular  residuals  as 
well.  Intelsat  4-F3  payload,  NORAD  id  5709,  yielded  one  of  the  poorest  fits.  Intelsat  4-F3 
payload  has  a  semi-major  axis  of  4.2349300e  +  04  km  and  a  period  of  8.6732156e  +  04 
seconds.  It  isn’t  in  the  close  vicinity  of  the  1:1  resonance.  Thus,  third  body  mass  pertur¬ 
bations  is  the  root  cause  of  the  poor  fit.  Figure  45  represents  the  position  and  the  velocity 
residuals  of  Intelsat  4-F3  payload  for  the  least  squares  fit.  The  propagation  of  Intelsat  4-F3 
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payload  from  01-01-2013  to  03-01-2013  yielded  an  rms  value  of  923.06  km. 
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Figure  44.  The  Residuals  for  Ops  9438  payload,  NORAD  id  10001,  over  2  Months 


4.4  Further  Analysis  for  Orbital  Characteristics 

The  low  eccentricity  KAM  theory  is  known  to  fail  in  the  close  vicinity  of  the  polar  and 
the  critical  inclinations  and  eccentricities  above  1 0  3 .  However,  there  is  no  experimental 
proof  for  the  limits  of  the  theory  in  terms  of  inclination  and  eccentricity.  This  section 
rigorously  investigates  those  limits.  For  this  section,  890  new  test  cases,  which  have  an 
eccentricity  of  1 0  3  and  smaller,  for  polar  inclination,  432  new  test  cases,  which  have  an 
eccentricity  of  10  3  and  smaller,  for  the  critical  inclination,  and  160  new  test  cases  for 
eccentricity  were  retrieved  from  www .  space-track .  com  servers,  and  propagated  from  01- 
01-2014  to  03-01-2014.  Another  approach  has  been  taken  for  the  pseudo-random  selection 
of  the  new  test  cases.  For  each  case,  the  interval  of  interest  was  divided  into  4  equal  pieces. 
Equal  number  of  random  test  cases  were  selected  for  each  interval. 
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Figure  45.  The  Residuals  for  Intelsat  4-F3  payload,  NORAD  id  5709,  over  2  Months 


Figure  46  represents  the  success  ratio  versus  the  eccentricity.  The  range  for  the  eccen¬ 
tricity  is  between  0.01  and  0.05.  This  plot  shows  the  upper  limit  of  the  eccentricity  that  can 
be  modeled  by  the  theory.  The  drastic  decrease  in  the  success  ratio  can  be  seen  in  the  plot. 
The  least  squares  fit  clearly  fails  for  the  eccentricities  of  0.0233  and  bigger.  The  decrease 
in  the  success  ratio  for  the  eccentricity  from  0  to  0.01  is  25%,  see  Section  4.2.  However, 
the  decrease  in  the  success  ratio  for  the  eccentricity  from  0.01  to  0.014  is  approximately 
45%. 

Figure  47  represents  the  success  ratio  versus  the  inclination,  which  is  in  the  close  vicin¬ 
ity  of  the  polar  inclination.  From  84°  to  96°  all  311  test  cases  have  failed  to  converge. 
Therefore,  the  close  vicinity  of  the  polar  inclination  certainly  terminates  the  fitting  process 
completely.  The  current  theory  can’t  model  orbits  which  have  inclinations  between  84° 
and  96°.  As  mentioned  previously,  the  solution  to  this  problem  is  known,  but  it  is  left  to 
future  studies.  This  work  is  specifically  beneficial  for  collision  avoidance  calculations  and 
formation  flight  applications. 
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Figure  48.  The  Optimal  Region  in  the  Close  Vicinity  of  Critical  Inclination 
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Figure  48  represents  the  success  ratio  versus  the  inclination,  which  is  in  the  close  vicin¬ 
ity  of  the  critical  inclination.  The  interval  between  64.5°  and  59.5°  is  where  the  success 
ratio  becomes  zero.  Thus,  the  optimal  region  is  the  interval  between  65.5°  and  58°.  Al¬ 
though  it  is  hard  to  define  the  starting  point  of  the  optimal  region  due  to  the  missing  points, 
the  upper  limit  is  definitely  64.5°. 

The  author  also  analyzed  the  13:1  resonance.  156  tests  cases  failed  between  6570  sec¬ 
onds  and  6700  seconds.  Delta  1  debris,  NORAD  id  6213,  Jason  debris,  NORAD  id  35414, 
and  Cosmos  1691,  NORAD  id  35414,  were  the  only  survivors  in  the  close  vicinity  of  the 
13:1  resonance.  Figure  49  represents  the  success  ratio  versus  the  period.  This  analysis 
provides  valuable  information  for  the  behavior  of  the  new  theory  near  resonance. 

4.5  Comparison  between  KAM  Torus  Method  and  SGP4 

The  least  squares  fitting  process  of  the  low  eccentricity  KAM  torus  to  the  SGP4  and 
TLE  data  requires  continuous  differential  corrections  of  B*  value,  basis  frequencies,  and 
the  the  global  variables.  The  author  has  written  additional  scripts  to  compare  the  new 
method  to  SGP4.  First,  a  low  KAM  torus  is  built,  and  then,  least  squares  is  fitted  to  SGP4 
and  TLE  data  over  two  months.  Theory  files,  which  are  used  to  build  the  torus,  and  the 
basis  frequencies  are  stored.  Then,  TLE  of  a  future  date,  which  is  two  months  later  than 
the  starting  date  of  the  propagation,  is  propagated  for  an  orbit  by  SGP4.  This  propagation 
is  assumed  to  be  the  truth.  Next,  both  the  low  eccentricity  KAM  torus  and  SGP4  are 
propagated  over  two  months  and  residuals  are  calculated.  The  same  process  is  repeated  for 
the  second  time  for  other  TLEs. 

Figure  50  represents  the  residuals  over  one  period  for  the  new  theory.  Figure  5 1  rep¬ 
resents  the  residuals  over  one  period  for  SGP4.  This  comparison  is  made  for  the  best 
least  squares  fitted  near-Earth  object,  which  has  an  rms  value  of  0.377006  km  over  two 
months.  This  test  case  is  Thorad  Delta  1  debris,  NORAD  id  8168.  It  has  an  inclination  of 
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Figure  51.  SGP4  prediction  over  2  Months  for  the  Best  Fit  Case 
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1.01334e  +  02°,  a  semi-major  axis  of  7.83 195e  +  03  km,  and  B*  value  of  8.21740^  —  05. 
The  new  theory  yields  approximately  five  times  more  accurate  predictions  than  SGP4  does. 
The  theory  files  for  this  comparison  was  created  from  01-01-2013  to  03-01-2013.  The  date 
for  the  truth  TLE  is  05-01-2013. 


The  KAM  Theory  Prediction  For  Ops  6630 


Figure  52.  The  Low  Eccentricity  KAM  Torus  Prediction  over  2  Months  for  a  Mean  Case 

Figure  52  represents  the  residuals  over  one  period  for  the  new  theory.  Figure  53  repre¬ 
sents  the  residuals  over  one  period  for  SGP4.  This  comparison  is  made  for  an  average  least 
squares  fit,  which  has  an  rms  value  of  5.0332  km  over  two  months.  This  test  case  is  Ops 
6630  debris,  NORAD  id  9323.  It  has  an  inclination  of  9.7031200e  +  01°,  a  semi-major 
axis  of  7.8056200e  +  03  km,  and  B  value  of  2.5971  OOOe  —  04.  The  new  theory  yields 
approximately  three  times  more  accurate  predictions  than  SGP4  does.  The  theory  files  for 
this  comparison  was  created  from  01-01-2013  to  03-01-2013.  The  date  for  the  truth  TFE 
is  05-01-2013. 

The  author  also  propagated  the  best  fit  case,  NORAD  id  8168,  from  03-01-2013  to  05- 
01-2013  with  approximately  two  days  intervals.  The  result  is  very  promising  because  there 
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is  a  linear  error  growth  for  one  month  with  a  small  slope.  Then,  the  error  growth  nearly 
doubles  between  30  to  60  days.  Figure  54  represents  the  residual  growth  versus  period. 
The  inaccuracies  in  SGP4  and  TLE  is  also  uncovered  by  Figure  54. 
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V.  Conclusions  and  Recommendations 


This  chapter  presents  the  conclusions  and  the  recommendations  for  future  studies.  The 
first  section  introduces  the  conclusions  obtained  from  chapter  4.  The  limitations  and  appli¬ 
cability  of  the  new  theory  is  presented.  For  the  second  section,  recommendations  for  future 
studies  are  given. 

5.1  Conclusions 

Although  the  title  for  this  effort  is  “KAM  Torus  Orbit  Prediction  using  SGP4  and  TLE”, 
it  has  yielded  more  than  only  the  possibility  of  orbit  prediction  combining  the  observational 
data  and  the  theory.  The  new  theory  has  proven  to  be  a  better  substitute  for  SGP4.  Its 
numerical  set  of  algorithms  provide  higher  accuracy  and  increase  computational  speed. 
The  theory  can  be  implemented  in  less  than  a  minute  on  a  home  computer.  The  new  theory 
has  binary  files  associated  with  parameters,  such  as  the  periodic  orbit,  the  geopotential  and 
air  drag  perturbations,  and  second-order  eccentricity  perturbation.  The  new  theory  contains 
two  parts.  First,  a  data  package  which  includes  Fourier  series.  Second,  a  set  of  seven  values, 
Q\.  Q2,  V2,  );3,  >’4,  V5,  y6,  and  B*  at  an  epoch  t0.  Therefore,  the  new  theory  seems  to  be  an 
alternate  way  to  compress  ephemerides. 

The  low  eccentricity  KAM  method  appears  to  be  promising  for  collision  avoidance 
calculations  because  it  yields  accurate  predictions.  The  linear  residual  growth  of  the  new 
theory  for  short  time  intervals  proves  its  applicability  to  collision  avoidance  calculations. 
Although  it  is  hard  to  make  a  precise  accuracy  analysis  without  raw  observational  data,  the 
error  growth  plot,  see  Section  4.5,  proves  that  the  low  eccentricity  KAM  torus  theory  may 
rival  numerical  methods  in  terms  of  accuracy.  In  addition,  the  new  theory  can  yield  far 
better  residuals  with  two  simple  modifications.  First,  the  relatively  inaccurate  TLEs  can  be 
eliminated  from  the  fitting  process,  see  Figure  54.  Second,  the  linear  correction  behavior 
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for  the  frequencies  during  the  fitting  process  can  be  projected  forward.  Moreover,  the  new 
theory  can  reduce  the  dependence  on  NORAD.  Because  most  ground  stations  depend  on 
NORAD  TLE,  they  must  wait  NORAD  to  publish  the  TLE.  However,  debris  and  most 
satellites  can  be  propagated  forward  in  time  using  TLE  history  with  high  accuracy  using 
KAM  torus  theory.  Formation  flight  is  another  possible  application  for  the  low  eccentricity 
KAM  torus  theory.  Since  the  accuracy  of  the  theory  increases  for  short  time  periods,  the 
new  method  can  be  an  onboard  collision  avoidance  algorithm  given  that  the  satellite  can 
receive  position  data. 

5.2  Future  Studies 

There  are  some  issues  related  to  the  new  theory.  As  presented  throughout  chapter  4, 
resonance,  high  eccentricity,  polar  inclination,  and  critical  inclination  issues  should  be  ad¬ 
dressed.  This  effort  provides  an  extensive  analysis  of  the  new  theory.  The  generalization 
of  the  theory  to  all  orbits  can  be  possible  after  these  issues  are  solved.  However,  this  effort 
can  be  used  as  future  reference  before  applying  theory  on  different  debris  and  satellites. 
The  solution  to  the  the  polar  inclination  problem  is  already  known.  The  problem  emerges 
due  to  the  forms  of  the  Legendre  polynomial  recursions  which  becomes  numerically  un¬ 
stable.  The  solution  to  this  problem  is  provided  in  Wiesel’s  text  [64].  The  atmospheric 
model  should  be  improved.  Russian  GOST  Atmosphere  is  a  promising  candidate  because 
it  requires  less  computation  time  and  it  is  relatively  accurate.  It  is  built  by  empirical  data 
obtained  from  Cosmos  satellites.  It  includes  solar  flux  and  geomagnetic  activity  effects 
[60].  The  application  of  the  theory  to  orbits  with  high  eccentricity  and  critical  inclina¬ 
tion  is  an  area  of  current  interest.  In  addition,  the  actual  accuracy  analysis  requires  raw 
observational  data,  because  SGP4  and  TLE  has  intrinsic  inaccuracies. 
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